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Abstract 


Many of the defects associated with solidification occur because of solute segregation in an 
alloy. In this thesis, we propose a model to track the levels of solute in a copper-tin alloy 
solidifying under natural convective forces. The thermal energy associated with a phase 
change is accounted for by an enthalpy dependent formulation of the specific heat. The 
discontinuity of solute at the interface is resolved by recasting the solute as the chemical 
activity, while treating the interface as an activity source. Through a two-dimensional finite 
element simulation, we are able to observe the flow pattern that develops in the fluid. The 
advection of solute due to buoyancy forces directly impacts the shape and progression of 
the interface. This model can be further expanded to three-dimensional shapes and many 


different materials. 
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Chapter 1 


Introduction 


1.1 Casting Review 


Casting of metallic alloys is an intensive process in terms of the time and resources consumed. 
The first step is to melt the alloy so that it can be poured into the appropriate shape. As 
the alloy begins to solidify, its components do not distribute homogeneously throughout 
the casting. The amounts of the different components at the solid-liquid boundary have a 
profound effect on the alloy that is produced. In particular, for peritectic solutions, such 
as iron-carbon or copper-tin, the concentration of solute determines which phases will form, 
and segregation may cause defects that negatively impact the quality of the finished product. 

Several factors influence how the solute disperses. Typically, the solute distributes de- 
pending on how well it can diffuse into the solvent. However, differences in the fluid density 
can produce internal currents, called convection patterns, that can advect the solute through 
the domain as well. Such variations in density can be caused by thermal fluctuations in the 
fluid or the quantities of each component if they have dissimilar individual densities. 

In order to prevent these defects from forming, it is important to track the presence 
of each metal during the solidification process. In this thesis we develop a model to find 
the temperature, velocities, and solute composition during the freezing process. Chapter 
2 derives the balance equations used to solve for these variables. Next, in Chapter 3 we 
expand upon how these equations are incorporated into a finite-element formulation. Lastly, 
in Chapter 4 we validate the model against several known solutions, and present results for 


a natural convection simulation. 


1.2 Previous Work 


The study of solidification physics began with the freezing of single materials. Such studies 


of ‘pure’ substances started with observations of ice formations from liquid water. Stefan 


oy : 
Gg" 
. _ « | 






















esivalhl wTEAO 


: quent oo ert wh ees 
ai” ou GCE > ees e eat 

_ Cae rapes 

ACP Te 

ain ole n> TT 

& apes % Srey 

af? aoe.) UG Tar6 OP 

view Goa 4 

: alviy pen 0 tip weal. 

re wtnet’ ri Gow er at CA 


7 erin Cuil toe oat 


- jews psi jee oe oie 


- Rar® GG) teqmne 06 
qWitiee wie aelinty aay 
a 7 
A rT item ii ae ver ‘ie i 
: 
uc «€ ; a | iA/ Oe \i'' wna Aj 
: ee eee 


ye Wotan eel) Gh anaes FO 
male ieeeinte = therein 


; 

salle olga at 
0 aha ee ee | 
an — a 


had previously worked on the flow of thermal energy across a fixed boundary. However, this 
new heat transfer problem involved a moving boundary between the phases [1]. His work 
led to what is now known as the Stefan condition. It is used to model the energy transfer 
when a phase change takes place. 

After the study of pure substances, further research was done on the freezing of solutions 
of two dissimilar materials. This was typically studied through the combination of different 
metals, resulting in a binary alloy. When unequal masses of the two components are used, 
the lesser is called the solute, and the greater is the solvent. For most alloys, the solute is 
considered to be homogeneously mixed in the solvent. Tiller et al. calculated how the solute 
distributes during one-dimensional solidification when the interface travels at constant speed 
and there is no convection in the liquid phase [2]. 

Smith et al. examined the solute distribution in the different regions of freezing, such as 
the initial transient, steady-state growth, and final transients [3]. Through the use of Laplace 
transforms, the equation for solute transport was solved for the varying regions. Additional 
investigations were conducted for sudden changes in the solidification rate and for multiple 
solutes. 

The next step was to incorporate the thermal effects together with the solute into a 
single model. Rubinstein found a similarity solution for the case where the alloy freezes at 
a specific, but unspecified, temperature in a semi-infinite domain [4]. His analysis has a 
boundary that moves according to the thermal balance. 

After the derivations of analytical solidification solutions, numerical modeling of the 
solidification began to develop. However, there are several obstacles that make this difficult. 
For example, most materials have material properties that vary between the phases. Several 
different approaches were taken to account for these variations. Ferreira et al. use different 
densities, specific heats, and thermal conductivities in the solid and liquid phases [5]. Their 
model defines the material to be liquid if it is above the liquidus temperature and solid if it 
is below the solidus temperature. Between these two temperatures, an average of the solid 
and liquid material properties is used. 

Another way to incorporate the differences in material properties is through a weighted 
average of the values based on the amount of each phase. Two ways to evaluate the amount 


of material present are the mass fraction and the volume fraction. They are defined as [6] 


M=g9.M,+g9eMe ; V=fsVs+ feVe 
GQotge=1 ; fot fe=l (1.1) 
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where M is the total mass of the system, g is the mass fraction, V is the total volume, 
and f is the volume fraction. The subscripts s and @ represent the solid and liquid phases, 


respectively. The phase fractions are related to each other such that 


Ge = fepe; PGs = fePs pisZ) 


where p is the density. In cases where the density remains constant throughout solidification, 
the mass fraction and volume fraction are equivalent. For a flat and stable interface a linear 


interpolation can then be used for any material property so that 


w= fsvs+ fewe (p = constant) (1.3) 


where w is any material parameter. 

Another difficulty in the modeling of solidification is how to account for the interfacial 
conditions. Since the solid-liquid boundary moves through the domain, the latent heat 
associated with the change in phase cannot be applied to a fixed boundary. One way to 
account for this is to track the enthalpy changes in both phases simultaneously. The enthalpy 
of the solid is simply the sensible heat, while the liquid enthalpy is the sum of sensible heat 
and latent heat. The governing equation for heat can then be expressed as the sum of heat 
transfer from both phases [7]. 

The latent heat can also be treated as a heat source in the energy balance. Voller 
and Swaminathan derived several models that use constant values for the density, specific 
heat, and thermal conductivity in both phases, with an additional term to account for the 
difference in enthalpy of the two phases. The release of the latent heat is regulated by the 
rate at which the solid volume fraction changes over time [8-10]. This method conserves 
enthalpy, but spreads the interface over a finite region. 

A more precise way to account for the latent heat is to incorporate it with the specific 
heat. Voller et al. used this method by designing a model that has a piecewise continuous 
formulation of the specific heat [11,12]. Instead of using the nominal values of the specific 
heat, an effective specific heat is defined as the local change in enthalpy for a change in 
temperature. These models use a function of the change in liquid volume fraction over time 
to calculate the appropriate specific heat values. This method is much more accurate for 
transient solutions and converges easily in implicit formulations. 

The most difficult part of the modeling the solidification process is the interfacial condi- 


tions for concentration. At the interface, there is a large discontinuity in the concentration 
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as solute diffuses from the liquid into the solid. One way to apply the interfacial condition 
is to only have the interface present at prescribed positions. A node-jumping scheme uses a 
variable time step size so that the interface always lies on a line of nodes at a time step [13]. 
However, this scheme can be unstable for cases where the time step size is large. Also, it is 
only applicable for one-dimensional solidification. 

Crowley and Ockendon solved the problem with discontinuity in a manner similar to the 
enthalpy method for the energy balance [14]. They recast the solute transport equations in 


terms of the chemical activity Ag/¢ 





Cs/e =< (1.4) 


with C/¢ as the compositions, a constant b, and the slopes m,/¢ from the phase diagram. The 
activity is continuous across the solid-liquid interface. This method also requires a source 
term for activity at the interface. The amount of source added is based upon the rate at 
which the liquid volume fraction changes over time. This method conserves solute, but it 
spreads the interface across a control volume. Palle and Dantzig showed that this method 
can be coupled with the energy balance equation in multiple dimensions [15]. 

Flow in the melt may also play a large role in the solidification development the shape of 
the interface and how it progress. Worster and Wettlaufer showed how different microstruc- 
tures can develop depending on the fluid flow near the interface [16]. They categorize the 


flow pattern by the solutal Rayleigh number, given by 


_ PIG BsACL* 


Ra 
UDe 


(a6) 
where g is the gravitational acceleration, (3, is the solutal expansion, AC is the characteristic 
solute difference, L is the cavity length, yz is the dynamic viscosity, and Dy is the solutal 
diffusivity. This serves as a ratio of solutal buoyancy force to viscous forces. 

Moussa et al. studied how the flow pattern near a boundary affects solidification. Their 
work shows that near a wall or interstitial fibers, the flow pattern slows [17]. This will result 
in lower concentration areas as the solute is advected away by stronger flows away from the 
wall. Additionally, the interface deforms with the fluid flow into a convex shape towards the 
fluid. 

Convective flow during solidification is caused by local variations in the liquid density 


due to differences in temperature and composition. Peppin et al. experimented with the 
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formation of convective flows of ammonium chloride. They found that for low solidification 
rates jets of rising fluid known as ‘chimneys’ can form at the interface [18]. For high solute 
concentrations, there is a larger convective pattern near the interface than in the nominal 
fluid. Melnikov et al. studied the freezing of isopropanal and found several internal convec- 
tion cells [19]. The number and size of these cells varied with the solidification rate and the 
temperature boundary conditions. 

Other experiments were performed by Liu and Trivedi to study the final microstructure 
of a lead-tin mixture for different amounts of buoyancy currents [20]. By varying the size of 
the ampoule, the amount of convection that develops could be changed. Smaller ampoules, 
for which convection was minimized, resulted in a banded structure. Larger samples had a 
more random microstructure development. 

Recent work done by Kohler combined solidification experiments with a finite element 
model. In his Ph.D. thesis, Kohler cast copper-tin alloys by directional solidification at 
very slow rates to ensure that the interface remained flat and stable throughout the cooling 
process. Thus, the effect of natural convection could be observed on the developing mi- 
crostructure [21]. He showed that the size of the sample, and thus the amount of convection, 
directly impacted the amount and type of phase that formed. 

Skudarnov et al. captured the flow pattern for aluminum chloride during solidification 
using particle image velocimetry (PIV) techniques [22]. They studied the convection currents 
in a rectangular cavity with cooled sides and thermally insulated top and bottom. For all 
cases, flow velocity is greatest at the beginning of solidification when there is a large domain 
of fluid. As more solid forms, the intensity of the flow pattern decreases. They also found 
that for higher concentration solutions, a larger number of convection cell patterns developed. 

These experiments and models all help to track the solute concentration during solidifi- 
cation. However, many of the models use an approximate formulas such as the Lever rule or 
Scheil equation, to find volume fraction as a function of temperature. It would be much more 
beneficial to have a model that evolves the interface based solely upon the local equilibrium 
conditions at the solid-liquid boundary. The goal of this thesis is to model the solute pattern 
of an alloy with this kind of interfacial condition in the presence of diffusive and convective 


effects. 
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Chapter 2 


Solidification Physics 


2.1 Frames of Reference 


The balance equations during solidification are used to compute the velocity, temperature, 
and concentration. In order to account for the flow of material, there are two frames of 
reference used to describe these variables [6, 23]. 

An Eulerian or laboratory frame of reference is fixed in space. A Lagrangian frame of 
reference is fixed upon a specific material point. As the material begins to change shape, the 
Lagrangian frame of reference changes with it. The Eulerian frame is convenient for defining 
boundary conditions. The Lagrangian frame allows for a better description of material 
interaction. 

In order to account for changes in the material, the time derivatives of the parameters 


must be accounted for in each frame of reference. The velocity is defined such that 


Ox 


= (2.1) 


v= 





xX 


where Z is the position vector in the Eulerian frame, ¢ is the time, and JU is the velocity 
vector. The subscript X denotes a fixed material point. The material derivative is then used 


to define the time derivative of a parameter at any particular point in the Eulerian frame by 


ee 


Dy _dv| , avr] _ dv 
salt. Gale 


[Seta ePaper 





“i (a. Vu) (2.2) 


with V as the gradient operator. In the following sections, all of the balance equations are 


written in the Eulerian frame. 
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2.2 Mass Transport 


The mass balance is given by [6,23] 


D . 
= +p(¥-a) =0 (2.3) 


The density can be expanded in a Taylor series in temperature and concentration of the form 


O 0 
p= pot (T — Te) = a 
Tes 








+ ...Higher Order Terms (2.4) 
To,Co 
where 7p and Cp are the reference temperature and concentration, respectively, defining the 


reference density, po. ‘The volumetric expansion coefficients for each of these variables are 
defined such that [21] 
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(2.5) 
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For small changes in temperature and concentration, an accurate approximation of density 
can be obtained by truncating the series after the first order terms. The simplifies Eq. (2.4) 


to 
p = po|1 — Gr (T — To) — Bc (C — Co)| (2.6) 


Substituting Eq. (2.6) into the mass balance gives the form 


Di po | Pr (1 — 1g) — Pa (C —.Ca) |} 
lye 
+po [1 — Gr (T — To) - Bo (C - Co) (¥-#) =0 (2.7) 


For the solidification problems studied in this paper, it is assumed that density variations 


due to temperature and concentration are small. This is a valid assumption so long as 
Br @& = To) <r ; Bo (C — Co) << (2.8) 
which gives the Boussinesq approximation [6, 16,23]. This reduces Eq. (2.7) to the form 


700 + 9) (7-0) = (2.9) 
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Since the derivative of a constant is zero, dividing Eq. (2.9) by po, gives the equation 
V-v=0 (2.10) 


which is referred to as the continuity equation. 


2.3 Momentum ‘Transport 


The change in momentum over time is given by the sum of external forces acting on a body. 


This can be described in each phase by the equation 


U = = & 
ae : ait 
PD, Vp+V-7T+ pb (217) 


where p is the pressure, T is the viscous stress tensor, and 6 is the body forces vector [6, 23]. 
The left side accounts for the transient changes in velocity and the advection of momentum 
through material flow. The terms on the right represent the pressure gradient, viscous forces, 
and body forces, respectively. The difficulty in analyzing this equation is that in natural 
convective flow, the density and pressure do not remain constant throughout the alloy. 
Most liquid metals behave as incompressible Newtonian fluids. The viscous stress tensor 


for a Newtonian fluid can be defined by 
cP ye D, peasy 


where D is the rate of deformation tensor [23]. Within each phase, the viscosity is constant, 


which allows the viscous term of Eq. (2.11) to be expressed as 
V-r=uvo (2.13) 
The modified pressure in a hydrostatic field can be defined by 
P= pot polgld (2.14) 


where d is the distance from the reference pressure, po, in the direction of gravity [6, 23]. 
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Combining Eqs. (2.13) and (2.14) into Eq. (2.11) gives the partial differential equation 


, Di 
po [1 — Br (T — To) = Oe Caran) Dt 
Se VE VEU bond Liege 0) & (2.15) 


The Boussinesq approximation can be invoked here to eliminate the effect of density change 
on the left hand side. This allows the governing differential equation for momentum to be 
expressed as 


f 


Dv > Bie. i, 
poap = —Vp + wV?8 + poBr (T — To) G+ poBc (C — Co) Gg (2.16) 


2.4 Energy ‘Transport 


The energy balance is given by 
ee (kV7) Pre D+ poQ (2517) 


where h is the enthalpy, k is the thermal conductivity, T is the temperature, and Q is a vol- 
umetric heat source [6,23]. The left hand side represents the transient flow of heat and the 
advection of heat through movement of the material. The first term on the right accounts 
for the thermal diffusion. The last two terms on the right side of the equation represent 
the heat generation by viscous dissipation and internal heat generation, respectively. Vis- 
cous dissipation is heat generated by internal mechanical forces due to deformation. For 
the solidification problems examined in this thesis, it is assumed that this contribution is 
negligible compared to the advective and diffusive terms. It is also assumed that there is no 
internal heat generation. The thermal conductivity is assumed to be constant throughout 
the freezing process in each phase. These assumptions simplify the energy equation to [10] 
Dh = 

Loa ar kW?T (2.18) 

As the alloy transforms from the liquid to the solid phase, additional energy is released. 
This energy, called the latent heat of fusion, must be conducted away from the interface 
in order for solidification to progress. This imposes the interfacial condition that the heat 


flowing into the interface along with the heat released by solidification must be balanced by 
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the heat flowing away from the interface. This is expressed in the form 


k (VT 





gis v7 it) = poL jt" ii (2.19) 
8 L 


with m as a unit normal vector, Ly as the latent heat of fusion, and v* is the interface 
velocity [6, 10, 23]. 
The interface condition can be combined with the energy balance in each phase to provide 


a governing equation for the entire alloy. The enthalpy is defined as 


df 
h = | adr; (2.20) 


To 


where c, is the specific heat capacity and 7p is a reference temperature, usually taken to be 
298K. For the problems we consider, the specific heat is assumed to be equal and constant 


within each phase. This allows the enthalpy to be expressed in the form 
h = cy (T — To) + febs (2.21) 


Using the chain rule, the derivatives of enthalpy can be expressed such that 


Oh Oh OT 

at aT ab ees 
=— —_VT Aes 

Vh= 2 (2.23) 


Within each phase, the derivative of enthalpy with respect to temperature is simply the 
specific heat. However, at the solid-liquid interface, there is a discontinuity in the enthalpy 
due to the latent heat release. In order to incorporate the interfacial condition, a new 


variable, the effective specific heat, is defined by 


c fe=1 
Oh 
oy te oh 0<fe<l (2.24) 
Cp We 


The derivative, 0h/OT is undefined at the interface. This problem is addressed in the 


numerical implementation chapter. 
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Thus the governing differential equation for energy can be expressed by the form 


oe Ale ~ 
pote ak ee (2.25) 


This form of the energy equation is much more convenient for solidification modeling for 
several reasons. Firstly, it allows for the tracking of the solid-liquid interface as it progresses 
through the domain. Also, it allows the intrinsic variable, temperature, to be continuous 


across the interface. Additionally, the enthalpy, is conserved throughout the freezing process. 


2.5 Solute Transport 


The solute balance in each phase is described by the equation 


= =V.(p¥c) +R (2.26) 


where R is the solute generation [6,23]. Once again, the left side of the equation represents 
the transient change in composition and the advection of solute through the movement of 
material. The first term on the right hand side describes the transport of solute by diffusion. 
The last term on the right represents the generation of solute by internal chemical reactions,, 


which we neglect in this work. This simplifies Eq. 2.26 to 


DC = ~ 

— =V- (pvc) 2.27 

Di (2.27) 
The solid-liquid interfacial condition for solute describes the balance of solute diffusion 

from each phase. The amount of solute diffusing into the interface and the solute rejected 

at the interface must be balanced by the solute diffused away from the interface. This is 


expressed in the form 


> 


DiVG\ st DV Gla = Coan (2.28) 
8 L 

where C7 and C* are the solute compositions at the interface for the liquid and solid phases, 

respectively [6,23]. Thus the diffusion of solute into and away from the interface is related 

to the thermodynamic equilibrium condition from the phase diagram. 


An alternative for can be obtained by replacing the composition with the chemical ac- 
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tivity, which is a continuous function, defined in each phase by [15, 24] 
Ap=C,; A, = ee (2.29) 


with ko as the segregation coefficient defined as 


m 
kyo = — (2.30) 

Ms 
where m, and m, are the liquidus and solidus slope. At the solid-liquid interface, the chemical 
activities of the two phases are equal. In order to relate the composition to the chemical 


activity, a new variable called the capacitance, Ccap is defined such that 


OC 





cap = = Bol 
C P OA ( ) 
From Eq. (2.29), the capacitance is constant in each phase as 
Ceap,t = 1 > Ccap,s = Ko (2.32) 
A law of mixtures can be used to define the capacitance at any point in the domain 
Ccap = faced 7 TeCeap.t = (1 = fe) Ccap,s + fACcane (2.33) 
The composition can then be recovered from Eq. 2.29 
i= 6 a (2.34) 
Substituting the chemical activity into the differential equation for solute yields 
D ca A — = 
(CeapA) why (Dv (Cap) (2.50) 
Dt 
Expanding this expression using the product rule yields 
|B iee® DA + 2 a | 
A + Cop = V (DeapVA + DAV ceap (2.36) 
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By use of the chain rule on Eq. 2.33 the derivatives of capacitance are given by 


Osea denny: 








= PS Hil 
Ot df, Ot eat) 
=, ie, 
VCcap = —— 2.38 
Eq. (2.37) can be computed directly so that the time derivative of capacitance is 
Orn Ofe 
a = (Geant = Coane) OL (2.39) 
Similarly, Eq. (2.38) reduces to 
items = (Crag e on cane) Vie (2.40) 


In the problems considered in this thesis, for a solidifying binary alloy, the interface is 
assumed to move slowly enough that the interface remains relatively flat and stable. This 
reduces the spatial derivative of the liquid fraction to zero within each phase. However, the 
spatial derivative will be infinite at the interface, which will be accounted for in the numerical 
implementation chapter. 


The differential equation for the chemical activity can now be expressed as 


df, DA 


Cea wa -Can.s A— + Cap 
(Ccap, — Ccap,s) Ae P Dt 


aa (DecapV A) (2.41) 


The first term on the left is brought over to the right side of the equation, where it resembles 
a source of activity. Thus the governing differential equation for the chemical activity is 


solved as 


DA 


__DA df 
SED 


=V- (Decap¥ A) — (Ceap,e — Ceap,s) ar 


(2.42) 


This form of the composition equation eliminates the discontinuity of solute at the interface 


while conserving solute in the domain. 
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Chapter 3 


Numerical Implementation 


3.1 Numerical Methods 


Computational modeling of solidification has several specific and unique properties that 
require particular capabilities from a finite element program not commonly found in com- 
mercial packages. Most of the requirements revolve around the issue of incorporating both 
solids and liquids into a single continuum. The difficulty arises due to potentially large 
local variations viscosity, diffusivity, and other material properties. Additionally, the pro- 
gram must allow user-defined equations for how these properties are calculated throughout 
a transient simulation. 

There are other factors that are beneficial to the solution strategy, but not necessarily 
required to solve the differential equations. For example, variability in element size allows for 
greater accuracy in areas of large gradients in the solution variables. Additionally, control 
over time step size can improve both the rate of convergence and accuracy of the solution 
while providing reasonable computation times. A graphical user interface helps to effectively 
build the mesh required for the simulation. Also, providing post-processing calculations on 
the solution allows for estimates of additional information. 

The program chosen to model the solidification process was the FIDAP?™ finite ele- 
ment program. This program, currently owned by the ANSYS!™ Corporation, fulfills all 
of the necessary requirements to solve solidification problems. The program uses Fortran 
programming in the implementation of its user supplied subroutines. This allows control 
over any internal calculations for material properties and boundary conditions. The limita- 
tions imposed upon these calculations because of the software are described in the following 


sections [25]. 
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3.2 Distinguishing Phases 


In order to determine the material properties that depend upon temperature, concentration, 
and liquid fraction, the first step is to determine the phase of each element. The solution 
variables, i.e. temperature, velocity, and activity, are calculated at the integration points. 
For a four-node quadrilateral element, the domain of each element is defined by the nodes 
at each of the corners. Therefore the solution variables at the nodes and their locations in 
space must also be determined. 


The solutions at the nodes and at the integration points are related by 
[N] sia ak et ‘ey 


where [N] is the matrix of shape functions evaluated at the integration points, jasper is 
the column vector of one of the solution variables at the nodes or the nodal position, and 
{S'™} is the column vector of one of the solution variables or position at the integration 
points [26]. For a four-node linear quadrilateral element with four integration points, [V] is 


square and non-singular. Thus the solutions at the nodes can be found by 
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Figure 3.1: Upper Corner of the Copper-Tin Phase Diagram 
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Figure 3.2: Copper-Tin Activity Phase Diagram 
{Sage = [N]~* fevaal (3.2) 


Once the nodal coordinate and the solution variables at the nodes have been calculated, 
it can be determined whether the element is completely solid, completely liquid, or contains 
the interface. The temperature is calculated at the integration points by FIDAP?™. Ad- 
ditionally, a corresponding temperature can be predicted from the activity and the phase 


diagrams in Figures 3.1 and 3.2 such that 
TUE mir ACC Ts (3.3) 


where 7’; is the freezing temperature of pure copper. The indicator variable, @ is then defined 


as 
i ie aa Tnode ee, dba = rode eck mA “e T; (3.4) 


If ¢ is greater than zero, this implies that this point lies above the liquidus line on the 
phase diagram. Thus the alloy at this node is in the liquid phase. Conversely, if @ is less 


than zero, the point lies below the solidus line, and the alloy at this node is solid. If ¢ is 
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equal to zero at the node, the interface lies directly on the node. For the element to be 
completely in the liquid phase, @ at each of the nodes must be greater than zero. Similarly, 
if @ at all four nodes of the element is less than zero, the element is solid. 

When ¢ is positive at some nodes in an element and negative at others, the interface lies 
within the element. In this case, a marching squares scheme is used to determine which of 
the sides of the element the interface passes through [27]. First, the nodes are numbered 
according to the schematic in Figure 3.3. The sum of nodes in the liquid phase or ‘on’ is 


found by 
4 
Tot =e 2 eC 0 (3.5) 
I 


where node is the number of the node. Each value of 0,, corresponds to a different case for 
which sides the interface passes through and where the solid-liquid boundary occurs. Figure 
3.4 shows all possible cases for these values. The two cases that are not implemented, case 


5 and case 10, are instances where two interfaces are present in the element. 








Node 4 Node 3 
Integration Point 2 Integration Point 4 
Integration Point 1 Integration Point 3 

Node 1 Node 2 








Figure 3.3: Finite element with FiIDAP?™ numbering sequence 
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Case 1 


Liquid Area = 
0.5*Phi1*Phi2 


Liquid Area = 


0.5*(dy-Phi1) 
*(dx-Phi2) 


Liquid Area = 
0.5*Phi2 
*(dy-Phi1) 





(g) Ton = 8 


Liquid Area = 
0.5*dx*((2*dy) 
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Figure 3.4: Control volume representation for different interface crossings 
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Figure 3.5: Control volume representation with nodal values 

















Since the interface is stable and planar, it is assumed that these situations are not possible 
within these simulations. Once the sides through which the interface of the element passes 


are determined, the length of the element sides is found by 


Iy2 = \f (%2 - £1) + (yo — yi) (3.6) 


where £2 and yo correspond to the nodal positions of the node farther away from node | and 
x, and y; are the nodal positions of the node closer to node 1. 

The next step is to find the point along the sides where ¢ is equal to zero. Both temper- 
ature and activity vary linearly along the side. Along the side, the point where the interface 


crosses is denoted as ¢p. This position can be found by the equation 


2 . o1 
Ly2 





oi + go = 0 (3a) 


where ¢, and @2 are the values of @ at the nodes at the ends of the side. Solving for ¢@9 and 


inserting the appropriate values for @, and @2 yields 


Tye (Ty = T\ oie mA) 


Ce eae ee ae 


(3.8) 
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Once the interface intercepts have been determined, the area of the liquid region can be 
determined by simple geometric calculations. The volume fraction of liquid in the element 


can then be calculated as 
fro= — (3.9) 


where ay is the area of the liquid region area and a,, is the area of the entire element found 
by 


1 
det = VV AL gL 34 — (Ly + 134 — Ly — 135)” (3.10) 


where the subscripts denote the node number [28]. 


3.3. User Defined Subroutines 


Several of the material properties depend upon the liquid fraction in the element and its time 
derivative. However, these material properties are calculated only at the integration points. 
In order to align these methods, the element is split up into four individual control volumes 
as shown in Figure 3.6. The same methodology used to determine the liquid fraction of the 
entire element can be used to determine the volume fraction within each control volume. 
In order to calculate the volume fraction of each control volume, the solutions and posi- 
tions are needed at the midpoints of each side as well as at the centroid of the element. Since 
the elements are four-node quadratic elements, the solutions vary linearly along the sides. 
The solutions and positions at each of these points can be found by linear interpretation to 


be 
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Figure 3.6: Finite element division into appropriate control volumes 


An additional change is that instead of Eq. (3.9), the liquid fraction is found by 


ij (3.12) 
aCV 


where the control volume area, acy is found by the inserting the lengths of the sides of the 
control volume into Eq. (3.10) 

The effective specific heat from the enthalpy method can then be calculated using Eq. (2.24). 
By using the temperature and liquid volume fraction at the previous and current time steps, 


the effective specific heat can be approximated in the presence of the interface as 


eft _ dh | Oh/dt | hi; - hj 
2 Sah ieee eT 





(3.13) 


where the subscript 7 denotes the integration point/control volume number and the super- 
scripts n and n — 1 correspond to the current and previous time steps, respectively [6]. The 
enthalpies are calculated directly using Eq. (2.21). 

The source term adds activity only at the integration points. Once the volume fraction 


of each control volume is found, the contribution of the source at each integration point is 
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calculated as 


n fel; “ fel 


R = — (Ccape — Ccap,e) Al” —*———4, (3.14) 
i" — ¢| 
Instead of Eq. (2.26), FIDAP?™ solves the species equation in the form 
DC 


Therefore, Eq. (3.14) must be multiplied by the density in order to provide the correct value 
for the source term. This gives the equation 

n —l 
n Sel; z= fil; 


R= —Po (Ceap,e = Ceap,s) Al; 


3.16 
a Ae = Ae ( ) 


Several of the variable material properties can be calculated directly from the volume 
fraction at the current time step [10]. A rule of mixtures is used to find the appropriate 
value between the two constant values from each phase. This methodology can be used to 


find the capacity, diffusivity, and viscosity by 


Ceap Ccap,s Ceap,e 
De ales Ls rer ey IBY, 
be Hs He 
Ceap,s Ceap,l 
= (ee ie  e), eapctese, «D7 (3.17) 
Hs He 


Since the diffusivity in the solid is typically 10~* times smaller than that in the liquid 
phase, D, is approximated by 


eae 1), (3.18) 


Additionally, as the liquid volume fraction approaches 0, the viscosity approaches infinity. 


We approximate this condition by 


[ts = 10*py (3.19) 
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Chapter 4 


Simulation Results 


4.1 Analytical Solutions 


Before using the previously discussed model to solve natural convection problems, a set 
of simpler problems with known solutions was tested. In this way, the volume fraction 
formulation of the activity-composition relationship was verified. The important points to 
confirm were the appropriate activity generation, while remaining continuous across the 
interface. Once the model was proven against known solutions, it could then be used to 
solve problems for which there are no analytical solutions. 

We use two solidification problems with analytical solutions to validate the model. Smith 
et al. found a solution for the one-dimensional solidification of a binary alloy where the 
interface moves at constant velocity [3]. The solution is obtained in the frame fixed on the 


interface by the coordinate transformation 
2=0—U t (4.1) 


where x is the absolute length from the edge of the domain, z is the distance from the 
interface, and v* is the prescribed interface velocity. The solution for the activity contains 
three distinct solidification regions: the initial transient, the steady-state, and the final 


transient. In the initial transient region, the activity in the liquid is given by 


a aoe hee 1 — ko ( =) F (5) Le (5) 
= exp { —— } erfc | ——— } — ~erfc | ——— 
i Jenn Guadve addy SABYI IN) 9 Pe 2/Dil 
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and in the solid phase by 


Ao 1 /vuta 
A, =—|1 f:{ — 
Do | Bi (; =| 


ko (1 : Ko) Pigie 2ko — ] UT 
2ko — 1 ee : : 
+ (2kp — 1) exp ( D, erfe 5 D, (4.3) 


After an initial length of approximately 3.D;/kov*, the exponential terms of Eq. (4.2) and 











(4.3) go to zero. This corresponds to the point where the diffusion of activity away from the 
interface matches the progression of the interface. Eq. (4.3) reduces to A, = Ao/kio activity 
in the solid. The activity in the liquid maintains a constant shape ahead of the interface, 


given by 








1 — ko te 
Abe Aye il = 4.4 
e | ae ko exp ( ~) (4.4) 


Once the diffusion tail ahead of the interface reaches the end of the domain, activity can 
no longer diffuse away from the interface into the liquid. The solid activity thus greatly 


increases during this final region of solidification, and is given by 


ie ss (2n +1) Lye) exp el (4.5) 


At A 
: [i (m + ko) De 





(ae 


where L is the length of the domain. An accurate approximation to Eq. (4.5) can be obtained 
by taking terms up to n = 3, as for greater values of n, the exponential term approaches 
ZeYO. 

The second solidification problem with an analytical solution was solved by Rubinstein 
et al. [4]. In this problem, the interface moves according to the thermal solution in a semi- 


infinite domain. This leads to an interface velocity that depends on time in the form 


Pl, (4.6) 
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where the constant ¥ is found as the solution to the transcendental equation 





Be — ; Serf (7) (: < aad exp (7°) erfe )) | 


1k ky? ) k 
-i1—-, 1—k ex erfc 
| =D) (1 — ko) yexp (= D, ( oa _ 


where 7), is the boundary condition temperature and 7p is the initial temperature. A so- 


_ meAg 
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(4.7) 





lution exists when 7; lies between the solidus and liquidus temperatures for the nominal 
composition. The Rubinstein solution does not have any initial or final transient regions. 


Instead, it uses the solution based on an alloy that solidifies at the temperature 


T* = T, + (To — Th) ex (7) (1 - pelea exp (77) erfc ©) (4.8) 


The activity at the interface can then be found from the interface temperature. The 


interfacial activity is then obtained from the phase diagram 


re 
loa 


AS (4.9) 


The activity in the liquid phase is then given by 





A*—A 
Ay = Ap + 70 erfe( 


erfec (7 —) 


while the activity in the solid phase is equal to the interfacial activity everywhere. Note that 


a) (4.10) 


this solution does not conserve solute. 


4.2 1-D Finite Difference Model 


A finite difference model was developed as a simple way to test the solidification model before 
applying it to more complicated schemes. For this thesis, the equations for activity were 
solved using MATLAB/™. This allowed the activity at the nodes to be obtained in matrix 
form by direct calculation of the solute transfer equation. Additionally, the finite difference 
method provides fine control over mesh size and the time marching scheme. 

The easiest way to implement the transient solution is through an explicit time-stepping 


scheme. There are two convergence criteria that need to met for such a scheme. The Courant 
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Table 4.1: Material Properties 
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Number, Co, is the ratio of the distance the interface progresses in one time step to the nodal 
spacing. In order to prevent the interface crossing multiple nodes in a single time step, the 
Courant number must satisfy the condition 


Db ih (4.11) 
Az 


Coane 
where At is the size of the time step and Az is the distance between adjacent nodes. The 
Fourier Number is defined as the ratio of diffusion distance in a time step to the nodal 
spacing. For an explicit scheme, the grid Fourier Number must be maintained such that 
At 
With these criteria in mind the domain was chosen to have a length of 10mm with a 
nodal spacing of 0.1mm. The first node was placed half of the nodal spacing from the left 
boundary, and the last node was placed half of the nodal spacing from the right boundary, 
so that the representative control volumes associated with these nodes would have the same 
size as the control volume for each of the intermediate nodes. A time step size of At = 0.1s 
satisfies both criteria in Eqs. (4.11) and (4.12) 


The first derivative is evaluated at the control volume boundaries by the equations 
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The diffusivity and capacity are evaluated in the same manner. The diffusive term, which 
involves a second derivative, is then calculated at the nodes by the derivative of the product 
of diffusivity, capacity, and the first derivative 

aA aA 
oO Tye OA rs [Dccap Se | right ¥ [Decap dr | test (4 14) 
Ox anit Az 
The temporal derivative is simply the difference in activity between time steps divided by 


the time step size. 


aA i A= “Alger. 


a A (4.15) 


No-flux activity boundary conditions are applied at the boundaries through the use of 
‘ghost nodes’. An extra node is placed outside of each boundary of the mesh. However the 
activity equations are not solved at these nodes. Instead, the activity value at these nodes 
is set to the value of the adjacent node. This ensures that the first derivative of activity at 
the first and last nodes is zero at the boundary. 

The solidification process is simulated by moving the solid-liquid interface through the 
domain. Thus, the position of the interface is known a priori as a function of time. The 
solid fraction can be calculated as the ratio of the distance the interface has progressed into 
the control volume to the nodal spacing. The liquid fraction can then be found since the 
sum of the liquid and solid fractions must be one. 

For the Smith problem, the interface velocity moves at constant speed. In this simulation, 
the interface velocity was chosen to be 0.0lmm/s. The interface position is simply the 
product of this velocity with the time. The total time for solidification is then the length 
of the domain divided by the interface velocity, or 1000s. The activity distribution was 
recorded every 10s of simulation time. 

Figure 4.1 shows the computed initial transient region of the simulation, where the ac- 
tivity begins to build up at the interface, along with the analytical solution. Figure 4.2 
shows the steady-state region where the activity at the interface remains constant. In both 
of these regions, the activity solution from the model agrees with the analytical solution. 
At the interface, the finite difference model over-predicts the analytical solution because the 
activity source term adds solute in this control volume, and it has not yet diffused out to 
the adjacent nodes. 


The final transient region is shown in Figure 4.3. There is only an analytical solution for 
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Figure 4.1: Activity distribution for the initial transient of the finite difference Smith prob- 
lem. 
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Figure 4.2: Activity distribution during steady-state of the finite difference Smith problem. 
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Figure 4.3: Activity distribution for the final transient of the finite difference Smith problem. 


the solid activity in this section. Solid activity increases as the activity can no longer diffuse 
ahead into the liquid and begins to collect at the boundary. The solid activity matches up 
with the analytical solution except at the last node. The solutions differ here because the 
representative control volume of the last node does not fully capture the activity generation 
and diffusion for the last amount of solid. 

In the Rubinstein problem, the interface position moves according to the thermal solution. 
For these problems, the interface position is a function of time such that 
ay uaa (4.16) 

pep 

This gives a high initial interface velocity that decreases over time. The total solidification 


time for the domain is then given by 


(CpL? 
A4ky? 





told (4.17) 


In the finite difference simulation, of this problem the thermal solution was prescribed 


such that the boundary temperature was 872.37°C, the solidus temperature, and the initial 
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temperature of 972.37°C,, the liquidus temperature. This gives a value of 0.0307 for ¥, 
resulting in a solidification time of 303.2s. 

The Rubinstein solution has no transient regions because it assumes solidification at at 
a constant temperature based on the boundary conditions. This eliminates any transient 
regions at the beginning and end of solidification. Since the model conserves solute, it is 
more accurate in these areas than the analytical solution. In the steady-state regions, the 
model gives equivalent activity distributions to the analytical solution in the liquid phase 
while the analytical solution will under predict the activity solid. As seen in Figure 4.4, the 
interfacial positions correspond between the model and the analytical solution throughout 


the domain. 
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Figure 4.4: Activity distribution of the Rubinstein problem. 


4.3 1-D Finite Element Model 


Now that the model has been tested by the finite difference method, it can now be used in 
a finite element application. In these simulations, the interfacial position is determined by a 
combination of the temperature and activity solutions. A one dimensional problem restricts 


the solution by preventing the flow of solute through advection. Additionally, this allows the 
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solutions to be compared against the same analytical solutions presented in the last section. 

The previously developed implicit scheme provides a convergent solution for most values 
of the Courant and Fourier Numbers. However, there are additional criteria necessary to 
ensure a stable and planar interface. As the interface moves through the domain, it develops 
an activity profile ahead of it as solute is rejected into the liquid. This liquid concentration 
can be converted to a liquidus temperature through the phase diagram. If the temperature 
gradient in the fluid ahead of the interface is lower than this liquidus curve, solid can form 
ahead of the interface [23]. To prevent this instability from occurring, the temperature 


gradient in the liquid must be chosen to satisfy the condition 


G> MeCp (1 — ko) v* 


4.1 


Additional convergence conditions must be met for the model to resolve the activity 
distribution in the liquid. As the activity distribution forms, a boundary layer develops 
ahead of the interface. The size of this boundary layer depends upon the diffusivity of 
activity in the liquid and the velocity of the interface. The mesh must be fine enough to be 


able to resolve this boundary layer. This gives the criterion of 


De 
<< 
Sie 


Ax 





(4.19) 


The finite element mesh was built with a nodal spacing of 0.2mm and a total length of 
5mm in the x-direction. In the y-direction, the mesh was designed to be only 2 elements 
wide. The boundary conditions could then be applied on the top, bottom, left, and right 
surfaces. This leaves a one-dimensional line of ‘free’ nodes through the center of the domain. 
The implicit time-stepping was designed with a maximum time step size of 0.1s. 

Another factor of convergence is how the interface moves between elements. Within a 
single time-step, the interface can move into and out of an element between iterations. This 
can destabilize the solution as the appropriate interface position cannot be determined. 

In order to provide a smoothly moving interface, a relaxation factor is used on the 


temperature and activity solutions of the form 
sit! = es 1 4+ (1-€)s8° (4.20) 


Typical values of € are 0.5 for temperature and 0.9 for activity. 


The interface was moved at a constant velocity through the boundary conditions on 
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temperature. Initially, the domain was given a constant temperature gradient. This tem- 
perature gradient was maintained during the entire simulation. The boundary temperatures 


were then lowered at a constant rate. The interface velocity is related to these conditions by 
v= [4ar20) 


where 7’ is the cooling rate. Keeping a constant temperature gradient and cooling rate 
ensures that the interface moves at a constant rate. No-flux activity conditions were applied 
to the boundaries. 

In order to capture all of the transient regions in the domain, the velocity was moved at 
a rate of 0.0lmm/s. From Eq. (4.18), the thermal gradient was chosen to be 250°C/mm, 
giving a cooling rate of 2.5°C/s. This gives a total solidification time of 500s. The activity 
distribution was again recorded every 10s of simulation time. 

The finite element solution agrees with the analytical solution everywhere except the 
element in which the interface is present. Figures 4.5 and 4.6 show this progression through 
the initial transient and steady-state regions, respectively, with the analytical solutions. 


Similar to the finite difference model, the activity does not agree in the element transferring 
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Figure 4.5: Activity distribution for the initial transient of the finite element Smith problem 
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Figure 4.6: Activity distribution during steady-state of the finite element Smith problem 
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Figure 4.7: Activity distribution for the final transient of the finite element Smith Problem. 
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from the liquid to solid phase as activity is generated within this element. Additionally, 
during the final transient region shown in Figure 4.7, the model reproduces the final transient 
activity in the solid, except for the last element to solidify. 

Another interesting analytical solution can be found by perturbing the steady-state region 
of the solution by a sudden increase in the interfacial velocity to another constant speed. At 
the solute cannot initially diffuse away from the interface fast enough. This causes another 
transient region where solute begins to build up at the interfacial boundary. As the interface 
continues to progress, the solutal diffusion rate will increase until another steady-state region 
develops. 

Smith et al. derived an analytical solution to this problem [3]. The amount of solute 
that builds up at the interface is a function of the ratio of the two interface velocities. ‘The 
solution is given as a function of the distance from the interfacial position where the change 


in velocity occurred. The activity in the solid is found by the equation 
A, = Ao c ~ 5 erfe RES 
+ (? ) i) oe (HS) enh AS x) 
+ p= j ices aa (Pe Rs) ay (4 =) 


where v3 is the second interface velocity, r is the ration of the initial interface velocity to 

















(4.22) 





the second interface velocity, and ¢ is the distance from the interface position at the velocity 
change. 

For this simulation, the initial interface velocity was set to 0.0lmm/s. Once the in- 
terface had reached the steady-state region, at 300s, the interface velocity was increased to 
0.02mm/s by doubling the cooling rate of the boundaries while keeping the thermal gradient 
constant. As shown in Figure 4.8, the shape of the activity solution matches up with the 
analytical solution. The numerical solution slightly under predicts the ‘bump’ in activity in 
the solid. But, it does recover to the new steady-state region in the appropriate amount of 


time. 
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Figure 4.8: Activity distribution for the change in interface speed problem. 


4.4 2-D Finite Element Model 


Two-dimensional solidification problems become much more complicated to solve than one- 
dimensional problems because the flow of liquid must be considered, which requires that 
we solve the momentum equation. This adds three new degrees of freedom: velocity in the 
x-direction, velocity in the y-direction, and the pressure. Additionally, the interface may no 
longer remain flat due to the advection of solute near the interface. 

In order to account for the pressure term in the momentum equations, a penalty method 
is used. This scheme does not solve the continuity equation. Instead, the pressure term is 


equated to 


t\es 
p=- (=) Vow (4.23) 


where 77 is the penalty parameter, set to 10~°. The pressure term is thus eliminated from 
the momentum equation and pressure is recovered for the velocity solution. For a four-node 
quadratic elements with four integration points, FIDAP?™ the pressure is computed at the 


area centroid of the element, allowing the pressure to be discontinuous between elements. 
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For the natural convection problem, the mesh needs to be much finer than those previ- 
ously used. In order to solve the momentum equation accurately near the boundaries, several 
nodes have must be within the boundary layer. Fifty nodes were place along a boundary 
of width 4mm in the horizontal direction. Such a small nodal spacing was used in order to 
capture the large changes in velocity in the x-direction. The vertical direction had 100 nodes 
on the 5mm height giving a nodal spacing of 0.05mm. This allows the non-planar interface 
to be approximated using linear elements. 

Similar boundary conditions were used in the two-dimensional problem as in the one- 
dimensional problem. A constant temperature gradient of 250°C/mm was applied with a 
constant cooling rate of 2.5°C/s to move the interface at an approximate rate of 0.01mm/s. 
No-flux boundary conditions were applied to all sides for activity. Additionally, no slip 
boundary conditions were applied to the walls by setting the velocity to zero in both the x- 
and y-directions. In order to control the velocities to obtain a converged solution, the mag- 
nitude of gravity was limited to 1% of its nominal value. Also, since the thermal boundary 
conditions were chosen to maintain interface stability, we consider only buoyancy forces due 
to solutal differences. 

The material closest to the vertical walls first begins to solidify since it is closest to the 
decreasing boundary temperature, and thus heat is extracted there first. Activity is added 
at this location according to the solidification mode. This action is analogous to solute being 
rejected into the liquid. Since the solute is less dense than the nominal fluid, it begins to 
rise, while the denser fluid sinks, generating a convection pattern in the liquid phase. 

Along the interface, the convection currents flow outward from the centerline towards 
the vertical walls. As the fluid flows in this direction, it carries the solute with it, increasing 
the activity near the walls, while decreasing it near the centerline. Increasing the activity 
lowers the temperature at which the fluid freezes. Thus, the fluid near the centerline freezes 
at a higher temperature than the fluid near the walls. This causes the interface to take on 
a more concave shape. Figure 4.9 shows how the interface position varies in the horizontal 
direction. 

The flow pattern in the fluid is related to the Rayleigh number, which is proportional 
to the cube of the flow length (Ra = p|g| G,ACL?/wD,). As the solid begins to take the 
place of fluid in the domain, the characteristic flow length begins to decrease. Therefore, 
as solidification progresses, the intensity of the flow pattern decays. ‘Thus, the highest 
magnitude of velocity appears during the initial solidification process, and decreases over 


time. Figure 4.10 illustrates this behavior. Typically, the maximum velocities are found 


37 












= - 
Peg He Oe 
_ Gap eomaean® af? Oslo oo ab 


ian GY an ebwod eV Mate ai; 
wails (ee ereal «dt at 


5 tno the ees os tee ee? S Me 


= fT « 1 ier of iru: r- en hi a als ai -_ jy am se ‘ 
> & ws par rage rye iets @ mes eae 
ass a” oer lon “7 @| 
? ' ’ all i 
acum +be97] “ 
_. netenteé 
be, « < Cuter 
ve vif v 7 ; 
y 
. ae art? 
in ae 
, : 
7 
50) ’ 7 1} 4 
a io ’ ' ye mye wy wi serena 7 
Oe Gal by y pe ay | ¢~ant awh oft ) 
yee | wy 6. Alii wl . : pat -_ (evel oa aud 1 


a 
ah hy ntl i aia Fae Gb) hs Timer) at! ae nikal | ry 


ah, wi > beenergtg Galen Dae Wis “) yr 1" “ere % aaa 
ia®, wap ag ote euwiieab Yn toe 2. ee : 





2 e F e—© Centerline 
a ~—« Midway 
aa Edge 


Distance From the Coldwall (mm) 





100 200 300 400 500 


Time (s) 


Figure 4.9: Interface position for different horizontal positions 
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Figure 4.10: Maximum velocity in the domain during natural convection. 
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at the centerline of the domain in the middle of the liquid, where the two convective cells 
converge and have fully developed. 

Figures 4.12-4.21 display the solute pattern throughout the domain during solidification. 
Also included in these figures are the interface positions and the streamlines in the fluid. 
As can be seen, the activity does not distribute evenly in either the vertical or horizontal 
direction. This is to be expected in the vertical direction, as the solute begins to build up 
at the interface. As has been seen, the convection pattern flows towards the walls from the 
centerline in the vicinity of the interface. This increases the activity closer to the walls. 


Figure 4.11 shows the activity profile for several vertical cross-sections of the domain. 
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Figure 4.11: Activity profile in the vertical direction for different horizontal positions at the 
end of solidification. 
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Figure 4.12: 
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(a) Contour plot of activity with the interface 


(a) Contour plot of activity with the interface 


(b) Vector plot of velocity, 


Natural Convection Simulation Results at 50s 
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(b) Vector plot of velocity. 


Figure 4.13: Natural Convection Simulation Results at 100s 
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(a) Contour plot of activity with the interface 
and streamlines. 


Figure 4.14: 
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(b) Vector plot of velocity. 


Natural Convection Simulation Results at 150s 
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(a) Contour plot of activity with the interface 
and streamlines. 


(b) Vector plot of velocity. 


Figure 4.15: Natural Convection Simulation Results at 200s 
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(a) Contour plot of activity with the interface 
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(b) Vector plot of velocity. 


Figure 4.16: Natural Convection Simulation Results at 250s 
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(a) Contour plot of activity with the interface 


and streamlines. 


(b) Vector plot of velocity. 


Figure 4.17: Natural Convection Simulation Results at 300s 
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(a) Contour plot of activity with the interface | (b) Vector plot of velocity. 
and streamlines. 























Figure 4.18: Natural Convection Simulation Results at 350s 
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(a) Contour plot of activity with the interface 
and streamlines. 





























(b) Vector plot of velocity. 


Figure 4.19: Natural Convection Simulation Results at 400s 
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(a) Contour plot of activity with the interface 
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Figure 4.20: Natural Convection Simulation Results at 450s 
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(a) Contour plot of activity with the interface 
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(b) Vector plot of velocity. 
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Figure 4.21: Natural Convection Simulation Results at 500s 
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(b) Vector plot of velocity. 
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Chapter 5 


Conclusions 


This thesis began with a derivation of the balance equations given the assumptions taken 
for solidification of a binary alloy. In Chapter 3, we discussed how these equations could 
be applied to a finite-element model. Specifically, it was discussed how to incorporate the 
solid-liquid interface into the calculations of the capacity, diffusivity, viscosity, specific heat, 
and activity source. Then, in Chapter 4, we validated the model against several analytical 
solutions using both a finite difference and a finite element method. Lastly, we analyzed the 
results of a two-dimensional solidification simulation with natural convection. 

It has been shown that natural convection can affect the interface during solidification. 
By changing the amount of solute near the interface, the flow patterns alter the local freezing 
temperature. This influences the shape and movement of the solid-liquid boundary. As the 
solid forms, the internal flow pattern changes to accommodate the smaller domain. The 
combination of all these effects result in an irregular solute pattern in the solid. 

The next step for this model would be to analyze situations with higher convective flows. 
This would increase the advection of solute in the domain and would further affect the shape 
and progression of the interface. Additionally, the model could be expanded to a three 
dimensional simulation. This would provide accurate information about solute pockets or 
other defects that would form for more complex shapes. However, this would require a 
considerable amount of computing power to handle the increased degrees of freedom and 
number of nodes. Eventually, this model could be used to predict phase transitions for solid 


state transformations. 
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Appendix A 
MatLab Code 


A.1l smith_finite_difference.m 


function smith_finite_difference(dx,L,dt,t_end, Vp) 
%, Create Mesh 

macy =)(0°dx:L)? 

x_node = ((.5*dx):dx: (L-(.5*dx)))?’; 

size_x = size(x_node); 

7, Create Time Steps 

ta="(O:dtit_end)”; 

size_t = size(t); 

%, Initial Conditions 


AO = 10; 
CO = 10; 
g0 = 1; 


A_present = AO*ones(size_x(1,1),1); 
A future = AO*ones(size_x(1,1),1); 
C_present = CO*ones(size_x(1,1),1); 
C_future = CO*ones(size_x(1,1),1); 
g_present = g0O*ones(size_x(1,1),1); 
g_future = gO0*ones(size_x(1,1),1); 
x_star = 0; 

cap = zeros(size_x(1,1),1); 
De=—zeros(size x(1,1),;1); 
dAdx_cv_left = zeros(size_x(1,1),1); 
dAdx_cv_right = zeros(size_x(1,1),1); 
becv left = zeros(size_x(1,1),1); 
D_cv_right = zeros(size_x(1,1),1); 
cap_cv_left = zeros(size_x(1,1),1); 
cap_cv_right = zeros(size_x(1,1),1); 
dgdt = zeros(size_x(1,1),1); 

tran = *zeros(size_x(1,1),1); 
dDcapdAdxdx = zeros(size_x(1,1),1); 
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addition = zeros(size_x(1,1),1); 

t_graph = 101; 

% Phase Diagram Constants (Cu-Sn) 

Tf = 1084.87; %K Pure Cu Freezing Temperature 
Tp = 798; “%K Peritectic Freezing Temperature 
Ca = 13.5; %wt% Peritectic Alpha Concentration 
Cb = 25.5; %wt/ Peritectic Beta Concentration 
ml = (Tp-Tf)./Cb; %K/wt% Liquidus Slope 

ms = (Tp-Tf)./Ca; %K/wt% Solidus Slope 

kO = ml./ms; % Segregation Coefficient 

Tl = (ml.*CO)+Tf£; %K Liquidus Temperature 

Ts = (ms.*CO)+T£; %K Solidus Temperature 

4 Material Constants (Cu-Sn) 

rho = 0.00775; %g/(mm*3) Density 

mu = 0.0004; %g*mm/s Kinematic Viscosity 

k = 0.35*0.1; “%W/(mm*k) Thermal Conductivity 


cp = 0.516; %J/(g*K) Specific Heat 

Lf = 236.1290323; %J/g Latent Heat of Fusion 
Dl = 0.0041; %(mm*2)/s Liquid Diffusivity 

Ds = D1/10000; %(mm*2)/s Solid Diffusivity 


cap_s = kO; % Solid Capacitance 

cap_l = 1; % Liquid Capacitance 

h 

7 Stability Check 

Co = Vp*dt/dx; 

if(Co >= 1) 
fprintf(’Courant Number is too large!\n’) 
return 

end 

Fo = Dl*dt/dx/dx; 

if(Fo >= .5) 
fprintf(’Fourier Number is too large!\n’) 
return 


figi = figure; 

hold on 

fig2 = figure; 

hold on 

foren = Pisizeit(i;1)-1 
/ Graphing 
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if (ns=="1) 


end 


figure (figl) 
plot (x_node,A_present, ’bx-’) 
figure(fig2) 
plot (x_node,C_present, ’bx-’) 
% Save Data too 
fidi = fopen(’smith_FD_A.dat’,’w’); 
bore = —4 17 sizerx( 151) 
If. =size.x(1,1)) 
fprintf(fid1,’%d %d \n’, x_node(i), A_present(i)); 
else 
fprintf(fid1,’%d %d \n&\n’, x_node(i), A_present(i)); 
end 
end 
fclose(fidl); 
y} 
fid2 = fopen(’smith_FD_C.dat’,’w’); 
for eies 6¢siz65x (1,1) 
ihieic=sizeux(1.1)) 
fprintf(fid2,’%d %d \n’, x_node(i), C_present(i)); 
else 
fprintf(fid2,’%d %d \n&\n’, x_node(i), C_present(i)); 
end 
end 
fclose(fid2) ; 
vf 
fid3 = fopen(’smith_FD_t.dat’,’w’); 
Pporintt (fidsyesd \neyet(ns1)); 
fclose(fid3) ; 


/ Left Shift 
A_present = A_future; 
C_present = C_future; 
g_present = g_future; 
4 Move Interface 
x_Star = Vp*t(nt1); 


Lor 


1 o= desize_x61F1) 

% Volume Fraction 

if(x_star <= x_cv(i)) 
g_future(i) = 1; 

elseif (x_star >= x_cv(it1)) 
g_future(i) = 0; 
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end 
for 


else 
petoture (i) i=01=(x_star-x_cv(i))/(eicvG+ti)-x_cv(i)): 
end 
4 Control Volume Constants 
/ Capacitance 
cap(i) = ((1-g_present (i))*cap_s)+(g_present (i) *cap_1l) ; 
4 Diffusivity 
D(i) = ((1-g_present(i))*Ds)+(g_present (i) *D1) ; 


i=l size rx (1,1) 
/, Spatial Gradients 
if (i==1) 
dAdx_cv_left(i) = 0; 
dAdx_cv_right(i) = (A_present (i+1)-A_present (i) )/dx; 
D_cv_left(i) = Ds; 
D_cv_right(i) = ((1-(.5*(g_present (i)+g_present (it+1))))*Ds) 
+((.5*(g_present(i)+g_present (i+1)))*D1) ; 
cap_cv_left(i) = cap_s; 
cap_cv_right(i) = ((1-(.5*(g_present (i)+g_present (it+1))))*cap_s) 
+((.5*(g_present (i)+g_present (i+1)))*cap_l) ; 
elseif (i==size_x(1,1)) 
dAdx_cv_left(i) = (A_present(i)-A_present (i-1))/dx; 
dAdx_cv_right(i) = 0; 
D_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present(i))))*Ds) 
+((.5*(g_present (i-1)+g_present(i)))*D1) ; 
D_cv_right(i) = Dl; 
cap_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i))))*cap_s) 
+((.5*(g_present (i-1)+g_present (i)))*cap_1); 
capleverignt (7), = rcapsly 
else 
dAdx_cv_left(i) = (A_present(i)-A_present (i-1))/dx; 
dAdx_cv_right(i) = (A_present(i+1)-A_present (i) )/dx; 
D_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i) )))*Ds) 
+((.5*(g_present (i-1)+g_present (i) ))*D1) ; 
D_cv_right(i) = ((1-(.5*(g_present (i)+g_present (it+1))))*Ds) 
+((.5*(g_present (i)+g_present (i+1)))*D1) ; 
cap_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i) )))*cap_s) 
+((.5*(g_present (i-1)+g_present (i)))*cap_1) ; 
cap_cv_right(i) = ((1-(.5*(g_present (i)+g_present (i+1))))*cap_s) 
+((.5*(g_present (i)+g_present (i+1)))*cap_1) ; 
end 
dDcapdAdxdx(i) = ((D_cv_right (i)*cap_cv_right (i) *dAdx_cv_right (i)) 
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end 


-(D_cv_left (i) *cap_cv_left (i) *dAdx_cv_left(i)))/dx; 


/ Temporal Gradients 

dgdt(i) = (g_future(i)-g_present(i))/dt; 

% Transient Term 

tran(i) = (cap_l-cap_s)*A_present (i) *dgdt (i) ; 

4 Future Values of A 

addition(i) = (dDcapdAdxdx(i)-tran(i))*dt/cap(i) ; 
A_future(i) = addition(i)+A_present (i) ; 


C_future(i) 


A_future(i)*(((1-g_future(i))*cap_s)+(g_future(i)*cap_1)) ; 


/ Graphing 

if(n == t_graph) 

t_graph = t_graph+100; 

figure (fig1) 

plot (x_node,A_present, ’bx-’) 

figure (fig2) 

plot (x_node,C_present, ’bx-’) 

; Save Data too 

fidi = fopen(’smith_FD_A.dat’,’a’); 


end 


for 


end 


1f=$1 ¥size-x(1;1) 
if (ig=sizerx ist) ) 

fprintf(fid1,’%d vd \n’, x_node(i), A_present(i)); 
else 

fprintf(fidi,’%d %d \n&\n’, x_node(i), A_present(i)); 
end 


fclose(fidi); 


h 


PrdJesetopen(? smith FDicrdat? fia’ ) ; 


for 


end 


ie=eiesize x(1),1) 
if(i~=size_x(1,1)) 

fprintf (fid2,’/d Ad \n’, x_node(i), C_present(i)); 
else 

fprintf(fid2,’%d 4d \n&\n’, x_node(i), C_present(i)); 
end 


fclose(fid2) ; 


i 


fid3 = fopen(’smith_FD_t.dat’,’a’); 
Pprintlerids ,wideineeetin;1)); 
fclose(fid3) ; 
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% Last Time Step 
if(m == size_t(1,1)-1) 


end 
end 


figure(fig1) 
plot (x_node,A_future, ’bx-’) 
figure (fig2) 
plot (x_node,C_future, ’bx-’) 
% Save Data too 
fidl*= fopen(’smith_FD_A.dat’,’a’); 
foreie=r.. size x(1,1) 
if (i*=size-x(1,1)) 
rprintfi (fidl, *4d 4d \n”) xinode(i), A_future(i)): 
else 
fprintettidien,d denen, x nodetl), Astuture ))). 
end 
end 
fclose(fid1); 
i 
fid2 = fopen(’smith_FD_C.dat’,’a’); 
for 2="i'sizeex(i ,1) 
‘ft (r=sizerx(is1)) 
Supine id2 7d 7.0. \0 .x_ node(1). |Ceruture 1)»; 
else 
poner ids .d.,d \ne\n’, x noda(i)) CetutureGd)); 
end 
end 
fclose(fid2) ; 
hh 
fid3 = fopen(’smith_FD_t.dat’,’a’); 
Pprmintttrida,* ,ds\n 4) t(ntl;1)); 
fclose(fid3) ; 


figure (fig1) 


xlabel (’ 
ylabel(’ 


Length (mm) ’) 
Activity (/-)*) 


title(’Activity Distribution’ ) 
figure (fig2) 
xlabel(’Length (mm) ’) 


ylabel(’ 


Concentration (wt%)’) 


title(’Concentration Distribution’ ) 
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A.2. rs_finite_difference.m 


function rs_finite_difference(dx,L,dt,t_end) 
% Create Mesh 

meow =2(0:dx:L)?; 

x_node = ((.5*dx):dx:(L-(.5*dx)))’; 

size_x = size(x_node); 

/ Create Time Steps 

Boe © O2dt : tuend).2; 

size_t = size(t); 

% Initial Conditions 


AO = 10; 
CO = 10; 
gO = 1; 


A_present = AO*ones(size_x(1,1),1); 

A_future = AO*ones(size_x(1,1),1); 

C_present = CO*ones(size_x(1,1),1); 

C future = CO*ones(size_x(1,1),1); 

g_present = gO*ones(size_x(1,1),1); 

g_future = g0*ones(size_x(1,1),1); 

x_star = 0; 

cap = zeros(size_x(1,1),1); 

D = zeros(size_x(1,1),1); 

dAdx_cv_left = zeros(size_x(1,1),1); 
dAdx_cv_right = zeros(size_x(1,1),1); 
D_cv_left = zeros(size_x(1,1),1); 

D_cv_right = zeros(size_x(1,1),1); 

Gapecvi leit = zeros(size_x(1,1),1); 
cap_cv_right = zeros(size_x(1,1),1); 

dgdt = zeros(size_x(1,1),1); 

tran = zeros(size_x(1,1),1); 

dDcapdAdxdx = zeros(size_x(1,1),1); 

addition = zeros(size_x(1,1),1); 

Eograph = 101; 

7 Phase Diagram Constants (Cu-Sn) 

Tf = 1084.87; %K Pure Cu Freezing Temperature 
Tp = 798; “%K Peritectic Freezing Temperature 
Ca = 13.5; %wt% Peritectic Alpha Concentration 
Cb = 25.5; %wt% Peritectic Beta Concentration 
ml = (Tp-Tf)./Cb; %K/wt% Liquidus Slope 

ms = (Tp-Tf)./Ca; /%K/wt/% Solidus Slope 

kO = ml./ms; % Segregation Coefficient 

Tl = (ml.*CO)+T£; %K Liquidus Temperature 
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Ts = (ms.*C0)+T£; “ZK Solidus Temperature 
% Problem Defined Paramaters 

T_inf = Tl; %C Inital Fluid Temperature 
4TO = Ts-10; %C Coldwall Temperature BC 
TO = Ts; 

/% Material Constants (Cu-Sn) 

rho = 0.00775; %g/(mm73) Density 

mu = 0.0004; %g*mm/s Kinematic Viscosity 
k = 0.35; “W/(mm*k) Thermal Conductivity 


cp = 0.516; %J/(g*K) Specific Heat 

Lf = 236.1290323; 4J/g Latent Heat of Fusion 
Dl = 0.0041; %(mm*2)/s Liquid Diffusivity 

Ds = D1/10000; %(mm*2)/s Solid Diffusivity 


cap_s = kO; % Solid Capacitance 
cap_l = 1; % Liquid Capacitance 
/ Dimensionless Groups 
alpha = k/(rhoxcp); %(mm*2)/s Thermal Diffusivity 
theta_f = (Tf-TO)./(T_inf-TO); % Dimensionless Equilibirium Freezing Temperature 
Ste = cp*(T_inf-TO) ./Lf; 
Lew = alpha./D1; 
/% Analytical Constants 
Phi_Equation = @(phi) (((theta_f-(erf (phi) .*(1+(sqrt(pi()) 
.*phi.*exp((phi.*2)).*erfc(phi) ./Ste)))).*(1-(sqrt (pi( .*Lew) .*(1-k0) 
.*phi.*exp(Lew.*(phi.*2)).*erfc(sqrt (Lew) .*phi))))-(m1.*CO./(TO-T_inf) )) ; 
phiO = fzero(Phi_Equation,0) ; 
hs 
% Stability Check 
Fo = Dlxdt/dx/dx; 
Petro 2=1%5) 
fprintf(’Fourier Number is too large! \n’) 
return 
end 
hh 
hh 
figi = figure; 
hold on 
fig2 = figure; 
hold on 
Pomn = bisize(t(1;1)-1 
/ Graphing 
if(n == 1) 
figure (fig1) 
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plot (x_node,A_present,’bx-’) 
figure (fig2) 
plot (x_node,C_present, ’bx-’) 
% Save Data too 
Pivie-eropent’ rs FD UA. dat’, °?w?); 
POL sie= 21 Csizerx(ist) 
Pieie—s) 2erx (11) ) 
Sprintt (fidi;’ 4d 4d \n*, x node(i)),yAspresent(i)); 
else 
fprintf(fid1,’%d %d \n&\n’, x_node(i), A_present(i)); 
end 
end 
fclose(fid1) ; 
y 
fid2.= tf open rs kDecedaté piw 
for te="trsizetx (1; t) 
if.(i7=sizesx(1,1)) 
forintt (fid29a/d (7d *inwisexenode (i) })Gepresent (i)) ; 
else 
fprintf(fid2,’%d %d \n&\n’, x_node(i), C_present(i)); 
end 
end 
fclose(fid2) ; 
he 
fid3 = fopen(’rs_FD_t.dat’,’w’); 
fprintf(fid3,’fd \n’, t(m,1)); 
fclose(fid3) ; 
end 
ZaLett. SHEL 
A_present = A_future; 
C_present = C_future; 
g_present = g_future; 
% Move Interface 
x_Star = 2.*phi0.*sqrt(alpha.*t(nt1)); 
for 4 = Pasize xGis1) 
% Volume Fraction 
if (xeatare<=ixecv(1)) 
g_future(i) = 1; 
elseif (x_star >= x_cv(it1)) 
g_future(i) = 0; 
else 
g_future(i) = 1-(x_star-x_cv(i))/(x_cv(it1)-x_cv(i)); 
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end 
/% Control Volume Constants 
7, Capacitance 
cap(i) = ((1-g_present(i))*cap_s)+(g_present (i) *cap_1) ; 
/ Diffusivity 
D(i) = ((1-g_present (i) )*Ds)+(g_present (i) *D1) ; 
end 
fori = 1:size2x(1,1) 
% Spatial Gradients 
if (i==1) 
dAdx_cv_left(i) = 0; 
dAdx_cv_right(i) = (A_present(it+1)-A_present (i) )/dx; 
D_cv_left(i) = Ds; 
D_cv_right(i) = ((1-(.5*(g_present (i)+g_present (it+1))))*Ds) 
+((.5*(g_present (i)+g_present (it+1)))*D1) ; 
cap_cv_left(i) = cap_s; 
cap_cv_right(i) = ((1-(.5*(g_present (i)+g_present (i+1))))*cap_s) 
+((.5*(g_present (i)+g_present (i+1)))*cap_1l) ; 
elseif (i==size_x(1,1)) 
dAdx_cv_left(i) = (A_present(i)-A_present (i-1))/dx; 
dAdx_cv_right(i) = 0; 
D_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i))))*Ds) 
+((.5*(g_present (i-1)+g_present(i)))*D1) ; 
Decveright (i) -= D1; 
cap_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i))))*cap_s) 
+((.5*(g_present (i-1)+g_present (i)))*cap_1) ; 
capleeveripht (i) = capll 
else 
dAdx_cv_left(i) = (A_present(i)-A_present (i-1))/dx; 
dAdx_cv_right(i) = (A_present(i+1)-A_present(i))/dx; 
D_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i) )))*Ds) 
+((.5*(g_present (i-1)+g_present (i) ))*D1) ; 
D_cv_right(i) = ((1-(.5*(g_present (i)+g_present (i+1))))*Ds) 
+((.5*(g_present (i)+g_present (it+1)))*D1) ; 
cap_cv_left(i) = ((1-(.5*(g_present (i-1)+g_present (i))))*cap_s) 
+((.5*(g_present (i-1)+g_present (i)))*cap_1) ; 
cap_cv_right(i) = ((1-(.5*(g_present (i)+g_present (i+1))))*cap_s) 
+((.5*(g_present (i)+g_present (i+1)))*cap_1) ; 
end 
dDcapdAdxdx(i) = ((D_cv_right (i)*cap_cv_right (i) *dAdx_cv_right (i)) 
-(D_cv_left(i)*cap_cv_left (i) *dAdx_cv_left(i)))/dx; 
7% Temporal Gradients 
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dgdt(i) = (g_future(i)-g_present(i))/dt; 


/, Transient Term 


tran(i) = (cap_l-cap_s)*A_present (i) *dgdt (i); 
/% Future Values of A 


end 
7uGLr 
if(n 


addition(i) = (dDcapdAdxdx(i)-tran(i))*dt/cap(i) ; 
A_future(i) 
C_future(i) 


addition(i)+A_present (i) ; 
A_future (i) *(((1-g_future(i))*cap_s)+(g_future(i)*cap_1)) ; 


aphing 
== t_graph) 


t_graph = t_graph+100; 


figure (fig1) 


plot (x_node,A_present, ’bx-’) 


figure (fig2) 


plot (x_node,C_present, ’bx-’ ) 
i Save Data too 


end 
ela 
a (ri 


fidl = fopen(’rs_FD_A.dat’,’a’); 
for -1s—41 seize. xGh. 1) 


if (i g=size *x (Ci41),) 
fprintf(fidi,’%d %d \n’, x_node(i), A_present(i)); 
else 
fprintf(fid1,’%d %d \n&\n’, x_node(i), A_present(i)); 
end 
end 


meLosectidi )- 
i 


bidze=etopen( rs_FD'C. date ya); 
foreue= -iasize-x(1,1) 
bEGie—sizerx (11) 
fprintf(fid2,’%d %d \n’, x_node(i), C_present(i)); 
else 
fprintf(fid2,’%d %d \n&\n’, x_node(i), C_present(i)); 
end 
end 
fclose(fid2) ; 


h 


fid3 = fopen(’rs_FD_t.dat’,’a’); 
Pind Ghide), 4 ,da\natotn, 1); 
fclose(fid3) ; 


st Time Step 
==) gize_t(1,1)-1) 
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figure (fig1) 
plot (x_node,A_future, ’bx-’) 
figure (fig2) 
plot (x_node,C_future, ’bx-’) 
%, Save Data too 
Piolemeropan( reekDYAvdat’ .a yi 
PuUreia=s | 5 izo x Cle) 
if (i7=size_x(1,1)) 
POrIncL (tid. /,0),da\nem xeno0de (1) AtuLure( 1). 
else 
forantt (Lidl. 2 /de/dexng\nw x node(i), A future (1) )>; 
end 
end 
fclose(fid1) ; 
yf 
fid2t=vfopen(’rs_FD.C. dat’ ,’a’); 
fon Ie=al i size-x( 171) 
if Gi" =size.x(1,1)) 
Porintie1d2e), ded e\i ee xenode (1) Cafiture (i)),) >; 
else 
PoLIner idl /d9,d0 \n&\n’. x node (i), sCeruture()) 
end 
end 
fclose(fid2) ; 
h 
Pidse—ptopen (rs.FDst.dat’,’a’); 
Borer (Pda ya,.ds\ne, t(ntl,1)); 
tcelosattid3s); 
end 
end 
figure (fig1) 
xlabel(’Length (mm) ’) 
weabel (Activity (-/-)’) 
title(’Activity Distribution’ ) 
figure (fig2) 
xlabel(’Length (mm) ’) 
ylabel(’Concentration (wt%)’) 
title(’Concentration Distribution’ ) 
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Appendix B 
FIDAP Input Files 


B.1 Smith Solution Input File 


/$$$$$$$$$$ DESCRIPTION $$$$$$$S$$ 
/ Intial-Steady_State-Final Problem 
/ Time-Varying Temperature BC 

/ No Flux Activity BC 

/ Interface Velocity = 0.01 mm/s 

/ 

PILE C. ) 

1-D DIRECTIONAL SOLIDIFICATION OF COPPER-TIN SOLUTION 
/ 

/$$$$$$$$$S$ FI-GEN $$$$$$S$$S 
FI-GEN(  ) 

[RK KKK POINTS RK KKK 
POINT (ADD , COORDINATE, X=0, Y=-2) 
POINT (ADD, COORDINATE, X=5, Y=-2) 
POINT (ADD, COORDINATE, X=0, Y=2) 
POINT (ADD , COORDINATE, x=5, Y=2) 

[eR RRA LINES RRR 

POINT (SELECT, ID=1) 

POINT (SELECT , ID=2) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=3) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT, ID=1) 

POINT (SELECT , ID=3) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=2) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

[eK KKKAKAK EDGING #****KKKKK 


62 





















AAS eoperret sects? 


paienes IT, 


=> 
— 


_ 


hae Daeg in 
> a verre .s de 
i» aavaine al 
? > ¢ ota 
w I Hat : : 


penangeses, UckF ie 
sa o¢hee ee zeae 
4 ~*, &G, Trees 
o¥ poe , STABLG 
tof Ard, OFERTR 
7 qe Maar 7 
wom CCE 
(sat. 

(ct. 
{ - » ¢ 
irae 


CURVE (SELECT, ID=1) 

CURVE(SELECT , ID=2) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=25 , RATIO=0 , 2RATIO=0 , PCENTER=0) 
CURVE (SELECT, ID=3) 

CURVE (SELECT, ID=4) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=2 , RATIO=0 , 2RATIO=0 , PCENTER=0) 
eR RRR MESH kK 

CURVE (SELECT , ID=1) 

CURVE (SELECT , ID=4) 

CURVE (SELECT , ID=2) 

CURVE (SELECT , ID=3) 

MFACE (WIRE, LABEL="cusn" , EDG1=1 , EDG2=1 , EDG3=1 , EDG4=1) 
MFACE(SELECT , ID=1) 

MFACE (MESH , MAP ,ENTITY="cusn" ) 

/ 

END( ) 

/ 

/$$$$$$$$$$ FI-BC $$$$$$HSH$HS$ 

/ 

F1-BC(. ») 

BGADD (SELECT , EDGE, ID=1) 

BGADD (SELECT , EDGE, ID=3) 

BGADD (ADD , EDGE, ENTITY="wall" , INCLUDE) 
BGADD (SELECT , EDGE, ID=2) 

BGADD (ADD , EDGE, ENTITY="right" , INCLUDE) 
BGADD (SELECT , EDGE, ID=4) 

BGADD (ADD, EDGE, ENTITY="left" , INCLDUE) 
END( ) 

/ 

/$$$$$$$$$$ FIPREP $$$$H$HS$$HS 

/ 

FIPREP( ) 

/****K*K*KKK MESH DIMENSIONS ********** 
/ Mesh Length: mm 


S12 =: 5 
/ Mesh Width: mm 
$H = 4 


/*x***x*KX PHASE DIAGRAM CONSTANTS ****##444 
/ Pure Copper Freezing Temperature: C 

$Tf = 1084.87 

/ Peritectic Temperature: C 

$Tp = 798 
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/ Peritectic Alpha Concentration: C 


SCae=013875 
/ Peritectic Beta Concentration: C 
$Cb = 25.5 


/ Liquidus Slope: C/wt% 

$ml = ($Tp-$T£)/$Cb 

/ Solidus Slope: C/wt7, 

$ms = ($Tp-$Tf£)/$Ca 

/ Segregation Coefficient: - 

$kO = $ml/$ms 

/**** KK KKK CAPACITY-ACTIVITY RELATIONSHIP *****## #4 
/ Solid Capacity: - 


$cap_s = $k0 
meuigquid Capacity: - 
$cap_l = 1 


/*****K*KA MATERIAL CONSTANTS ** kX 
/ Density: g/(mm*3) 

$rho = 7.75E-3 

/ Thermal Conductivity: W/(mm*K) 

$k = 0.35 

/ Specific Heat : J/(g*K) 

$cp = 0.516 

/ Latent Heat of Fusion: J/g 

$Lf = 236.1290323 

/ Liquid Diffusivity: (mm*2)/s 


$D1l = 4.1E-3 
/ Solid Diffusivity: (mm*2)/s 
$Ds = 4.1E-7 


/xxxx «eK INITIAL ACTIVITY CONDITIONS *4*k#K KH 
Petnicial Concentration: wt/ 


$C_tO = 10 
geimitial*Activity: - 
$A_tO = 10 


/******K*KKK DERIVED QUANTITIES ********** 

/ Liquidus Temperature: C 

$T1 = ($m1*$C_t0)+$Tf 

/ Solidus Temperature: C 

$Ts = ($ms*$C_t0O)+$TFf 

/********** TNITIAL TEMPERATURE CONDITIONS **+#+#*** 4+ 
/Initial Temperature of the Left Wall: C 

$T_x0_tO = $Tl 

/Initial Temperature of the Right Wall: C 
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$T_xL_tO = $T1+1250 

/Temperature Gradient: C/mm 

$dTdx = ($T_xL_t0-$T_x0_t0)/$L 

/*******KKKX PROBLEM DEFINED CONSTANTS #2 
/ Interface Velocity: mm/s 


$Vp = 0.01 

/ 

/** AK PROBLEM DEFINITION *****k**x 
/ 


PROBLEM (2-D , TRANSIENT , NONLINEAR , NOMOMENTUM , ENERGY , SPECIES , SPECIES=1) 
EXECUTION (NEWJOB) 
PRINTOUT (NONE) 
DATAPRINT (NONE) 
/ 
/**k** KKK CONTINUEM ENTITIES *********x 
/ 
ENTITY (FLUID, NAME="cusn" , PROPERTY="cusn") 
7 
/*******KAKA BOUNDARY ENTITIES **#*** x 
/ 
ENTITY (PLOT , NAME="wall") 
ENTITY (PLOT , NAME="left") 
ENTITY (PLOT , NAME="right") 
/ 
/**k***KKKKK* SOLUTION STRATEGY *********X 
/ 
SOLUTION(S.S.=20, VELCONV=0.001,RESCONV=0.001) 
RELAXATION () 
eee tO. 07 5>.0, 0, 0, 0.95 
TIME (BACK , DT=1E-6 , NSTEPS=1E5 , VARIABLE=1E-3 , INCMAX=1.5, DTMAX=0.1,TEND=300) 
/ 
/***k**K*KKKA MATERIAL PROPERTIES ********4* 
Z 
DENSITY (SET="cusn" , CONSTANT=$rho) 
CONDUCTIVITY (SET="cusn" , CONSTANT=$k) 
SPECIFICHEAT (SET="cusn" , SUBROUTINE=7 ) 
BaexOetO, $1_xL_tO, $ml, $Tf, Scp, SLL, $L 
DIFFUSIVITY (SET="cusn" , SUBROUTINE=6 ) 
$ml, $Tf, $Ds, $D1, $cap_s, $cap_l 
CAPACITY (SET="cusn" , SUBROUTINE=4) 
$ml, $Tf£, $cap_s, $cap_l 
i 
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/** KK KAA SOURCE CONDITIONS *%* «x 

if 

SOURCE (SPECIES=1 , ENTITY="cusn" , SUBROUTINE=6) 

Sml, Sif, $cap.s, $cap_1, $rho, $A_to 

/ 

/*****KAKKK FLUID INITIAL CONDITIONS *****xe*ex 

/ 

ICNODE (TEMPERATURE, POLYNOMIAL=1 , ENTITY="cusn") 

Gtex0 tO, $dTdx,1,0,0 
ICNODE(SPECIES=1 , CONSTANT=$A_t0, ENTITY="cusn") 

/ 

/#* «44 TOP/BOTTOM WALL BOUNDARY CONDITIONS *##44« KK 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="wall") 

/ Activity Flux = 0 

/ 

/****K*K*KX RIGHT WALL BOUNDARY CONDITIONS *****#**** 
‘4 

BCNODE (TEMPERATURE , FSUBROUTINE, ENTITY="right") 

/ Activity Flux = 0 

/ 

/x* eee LEFT WALL BOUNDARY CONDITIONS #** #44 KX 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="left") 

/ Activity Flux = 0 


/ 

RENUMBER (PROFILE) 

END( ) 

/ 

/$$$$$$$$$$ FISOLV $$$$$$$$$S 
/ 

CREATE (FISOLV) 

/ 
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B.2 Smith Solution Restart File 


/$$$$$$$$$$ DESCRIPTION $$$$$$$$$$ 
/ Intial-Steady_State-Final Problem 
/ Time-Varying Temperature BC 

/ No Flux Activity BC 

/ Interface Velocity = 0.01 mm/s 

/ 

TETLE( ) 

1-D DIRECTIONAL SOLIDIFICATION OF COPPER-TIN SOLUTION 
/ 

/$$$$$$$$$$ FI-GEN $$$$S$HS$$ 
FI-GEN(  ) 

[eK KKKKAKK POINTS RK KKK KK 
POINT (ADD , COORDINATE, X=0, Y=-2) 
POINT (ADD , COORDINATE, X=5 , Y=-2) 
POINT (ADD , COORDINATE, X=0, Y=2) 
POINT (ADD , COORDINATE, x=5, Y=2) 
[RRR KK LINES RRR 

POINT (SELECT , ID=1) 

POINT (SELECT , ID=2) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=3) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=1) 

POINT (SELECT , ID=3) 
CURVE (ADD , SEGMENT) 

POINT (SELECT , ID=2) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

[eR KKK EDGING 28K KK 
CURVE (SELECT , ID=1) 

CURVE (SELECT , ID=2) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=25 , RATIO=0 , 2RATIO=0 , PCENTER=0 ) 
CURVE(SELECT , ID=3) 

CURVE (SELECT , ID=4) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=2 , RATIO=0 , 2RATIO=0 , PCENTER=0) 
[RK RK KA MESH RRR RK KK 

CURVE (SELECT , ID=1) 

CURVE (SELECT , ID=4) 

CURVE (SELECT , ID=2) 

CURVE (SELECT , ID=3) 
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MFACE (WIRE, LABEL="cusn" , EDG1=1 , EDG2=1 , EDG3=1 , EDG4=1) 
MFACE (SELECT , ID=1) 

MFACE (MESH , MAP ,ENTITY="cusn") 

/ 

END(  ) 

/ 

/$$$$$$$$$S$ FI-BC $$$$$$SS$$ 

/ 

FI-BC( ) 

BGADD (SELECT , EDGE, ID=1) 

BGADD (SELECT , EDGE, ID=3) 

BGADD (ADD , EDGE ,ENTITY="wall" , INCLUDE) 
BGADD (SELECT , EDGE, ID=2) 

BGADD (ADD , EDGE, ENTITY="right" , INCLUDE) 
BGADD (SELECT , EDGE, ID=4) 

BGADD (ADD, EDGE ,ENTITY="left" , INCLDUE) 
END() 

f 

/$$$$$$$$$$ FIPREP $$$$$S$$$$ 

pe 

FIPREP( ) 

/***kKKKKKA MESH DIMENSIONS ** kee 
/ Mesh Length: mm 


SLe="5 
/ Mesh Width: mm 
$H = 4 


/********** PHASE DIAGRAM CONSTANTS *#+#***4*** 
/ Pure Copper Freezing Temperature: C 

$Tf = 1084.87 

/ Peritectic Temperature: C 

$Tp = 798 

/ Peritectic Alpha Concentration: C 

$Ca = 13.5 

/ Peritectic Beta Concentration: C 

$Cb = 25.5 

/ Liquidus Slope: C/wt% 

$ml = ($Tp-$T£)/$Cb 

/ Solidus Slope: C/wt% 

$ms = ($Tp-$Tf£)/$Ca 

/ Segregation Coefficient: - 

$kO = $ml/$ms 

[x4 KKK CAPACITY-ACTIVITY RELATIONSHIP RH RKK 
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/ Solid Capacity: - 


$cap_s = $k0 
/ Liquid Capacity: - 
$cap_l = 1 


/********** MATERIAL CONSTANTS *** «44 
/ Density: g/(mm*3) 

$rho = 7.75E-3 

/ Thermal Conductivity: W/(mm*K) 

$k = 0.35 

/ Specific Heat : J/(g*K) 

$cp = 0.516 

/ Latent Heat of Fusion: J/g 

$Lf = 236.1290323 

/ Liquid Diffusivity: (mm*2)/s 


$D1 = 4.1E-3 
/ Solid Diffusivity: (mm*2)/s 
$Ds = 4.1E-7 


/****KAKKKA TNITIAL ACTIVITY CONDITIONS ********x 
Perhicial Concentration: wt% 


$C_tO = 10 
Peinitial “Activity: - 
$A_tO = 10 


/**k*k*KKKAKK DERIVED QUANTITIES ****+*+*4** 

/ Liquidus Temperature: C 

$Tl = ($m1*$C_t0)+$Tf 

/ Solidus Temperature: C 

$Ts = ($ms*$C_t0O)+$Tf 

/********** TNITIAL TEMPERATURE CONDITIONS ******k**x 
/Initial Temperature of the Left Wall: C 
$T_x0_tO = $Tl 

/Initial Temperature of the Right Wall: C 
$T_xL_tO = $T1+1250 

/Temperature Gradient: C/mm 

$dTdx = ($T_xL_t0-$T_x0_t0)/$L 

/********K** PROBLEM DEFINED CONSTANTS ***+*#&k**x 
/ Interface Velocity: mm/s 


$Vp = 0.01 

/ 

/x**kk*KKKK PROBLEM DEFINITION *##k44 44 
/ 


PROBLEM(2-D , TRANSIENT , NONLINEAR , NOMOMENTUM, ENERGY , SPECIES , SPECIES=1) 
EXECUTION (RESTART) 


69 












re es a 
| iw a 
odie) 
- We - 
1 sri Gere: 
we ecageet 
TERA J 
sree, 
3. Gee GATE q 
a” GAS be Ones 


. 
> oo ; 

morte & 

as en 4) iT. 

SOs 9 TS 

memene £79.) Grae re ia] ~ 


at : Ss 


-_ 


i ee te. | 





ery ages) rally Aare, Jai 


Wk 


a 


PRINTOUT (NONE) 
DATAPRINT (NONE) 
/ 
/***** KKK CONTINUEM ENTITIES ********** 
/ 
ENTITY (FLUID, NAME="cusn" , PROPERTY="cusn" ) 
/ 
/****KKKKAKA BOUNDARY ENTITIES ***k** kX 
/ 
ENTITY (PLOT , NAME="wall") 
ENTITY (PLOT , NAME="left") 
ENTITY (PLOT , NAME="right") 
/ 
/****K*KKKKK* SOLUTION STRATEGY ********KX 
/ 
SOLUTION (S.S.=20, VELCONV=0.001,RESCONV=0.001) 
RELAXATION () 
eeu). Oro, 0, 0,05, 0.95 
TIME (BACK , DT=1E-6 , NSTEPS=1E5 , VARIABLE=1E-3 , INCMAX=1.5,DTMAX=0.1,TEND=600) 
‘if 
/**k***K*KAKKA MATERIAL PROPERTIES ********** 
/ 
DENSITY (SET="cusn" , CONSTANT=$rho) 
CONDUCTIVITY (SET="cusn" , CONSTANT=$k) 
SPECIFICHEAT (SET="cusn" , SUBROUT INE=7 ) 
mrex0-tO, $T_xL_t0O, $ml, $Tf, $cp, $SLf, $L 
DIFFUSIVITY (SET="cusn" , SUBROUTINE=6) 
$ml, $T£, $Ds, $D1, $cap_s, $cap_l 
CAPACITY (SET="cusn" , SUBROUTINE=4) 
$ml, $Tf, $cap_s, $cap_l 
/ 
/[ **k**KKAKKKA SOURCE CONDITIONS ******4** 
/ 
SOURCE (SPECIES=1 , ENTITY="cusn" , SUBROUTINE=6 ) 
$ml, $Tf£, $cap_s, $cap_l, $rho, $A_tO 
/ 
/*********KX FLUID INITIAL CONDITIONS **** 444% 
/ 
ICNODE (TEMPERATURE, POLYNOMIAL=1 , ENTITY="cusn" ) 
miexOstoO, $dldx,1,0,0 
ICNODE(SPECIES=1 , CONSTANT=$A_t0O, ENTITY="cusn" ) 
if 


70 





areca G57 733 & Tree 
te Me, teec "Sa 


/ 444444 TOP/BOTTOM WALL BOUNDARY CONDITIONS **#4k#4«# 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="wall") 

/ Activity Flux = 0 

/ 

/ «HK RIGHT WALL BOUNDARY CONDITIONS #44444 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="right") 

/ Activity Flux = 0 

/ 

/*****K*KKK LEFT WALL BOUNDARY CONDITIONS *******xx 
, 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="left") 

/ Activity Flux = 0 


/ 

RENUMBER (PROFILE) 

END() 

“ 

/$$$P$$$$$$$ FISOLV $$$H$H$$S$SH$ 
/ 

CREATE (FISOLV) 

y 
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B.3 Change in Interface Velocity Input File 


/$$$$$$$$$$ DESCRIPTION $$$$$H$$$$ 
/ Intial-Steady_State-Final Problem 
/ Time-Varying Temperature BC 

/ No Flux Activity BC 

/ Interface Velocity = 0.01 mm/s 

/ 

BIILE( ) 

1-D DIRECTIONAL SOLIDIFICATION OF COPPER-TIN SOLUTION 
/ 

/$$$$$$$$$$ FI-GEN $$$S$$$S$$H$ 
FI-GEN(_) 

[RK KKK POINTS RK KKK 
POINT (ADD , COORDINATE, X=0, Y=-2) 
POINT (ADD , COORDINATE, X=30 , Y=-2) 
POINT (ADD , COORDINATE, X=0, Y=2) 
POINT (ADD , COORDINATE, x=30 , Y=2) 
[RR KKK LINES RRR KK 

POINT (SELECT, ID=1) 

POINT (SELECT , ID=2) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=3) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT, ID=1) 

POINT (SELECT, ID=3) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=2) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

[**kKKKKKAKKK EDGING *****KKKHX 
CURVE(SELECT , ID=1) 

CURVE (SELECT , ID=2) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=150, RATIO=0 , 2RATIO=0 , PCENTER=0 ) 
CURVE (SELECT , ID=3) 
CURVE(SELECT , ID=4) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=2 , RATIO=0 , 2RATIO=0 , PCENTER=0) 
[RK KKK MESH KKK KK 
CURVE(SELECT, ID=1) 
CURVE(SELECT , ID=4) 

CURVE (SELECT , ID=2) 

CURVE (SELECT , ID=3) 
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MFACE (WIRE, LABEL="cusn" ,EDG1=1 , EDG2=1 , EDG3=1 , EDG4=1) 
MFACE(SELECT , ID=1) 

MFACE (MESH , MAP ,ENTITY="cusn") 

/ 

END( ) 

/ 

/$$$$$$$$$$ FI-BC $$H$$H$$$HS$S 

/ 

FI-BC( ) 

BGADD (SELECT , EDGE, ID=1) 

BGADD (SELECT , EDGE, ID=3) 

BGADD (ADD, EDGE, ENTITY="wall" , INCLUDE) 
BGADD (SELECT , EDGE, ID=2) 

BGADD (ADD , EDGE, ENTITY="right" , INCLUDE) 
BGADD (SELECT , EDGE, ID=4) 

BGADD (ADD, EDGE, ENTITY="left" , INCLDUE) 
END(_) 

/ 

/$$$$$$$$$H FIPREP $$$$$SHS$HSH 

/ 

FIPREP(  ) 

/[*****kKKKKK MESH DIMENSIONS ****** KK 
/ Mesh Length: mm 


$L = 30 
/ Mesh Width: mm 
Shes) 4. 


/*********K PHASE DIAGRAM CONSTANTS *#****# 
/ Pure Copper Freezing Temperature: C 

$Tf = 1084.87 

/ Peritectic Temperature: C 

$Tp = 798 

/ Peritectic Alpha Concentration: C 

$Ca = 13.5 

/ Peritectic Beta Concentration: C 

$Cb = 25.5 

/ Liquidus Slope: C/wt% 

$ml = ($Tp-$Tf£)/$Cb 

/ Solidus Slope: C/wt% 

$ms = ($Tp-$Tf£)/$Ca 

/ Segregation Coefficient: - 

$kO = $m1/$ms 

/ 4K KKKK CAPACITY-ACTIVITY RELATIONSHIP 1K 
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/ Solid Capacity: - 


$cap_s = $k0 
/ Liquid Capacity: - 
$cap_l = 1 


/***** KA MATERIAL CONSTANTS *# #44 4k 
/ Density: g/(mm*3) 

Brno =.7.75E-3 

/ Thermal Conductivity: W/(mm*k) 

$k = 0.35 

/ Specific Heat : J/(g*k) 

$cp = 0.516 

/ Latent Heat of Fusion: J/g 

$Lf = 236.1290323 

/ Liquid Diffusivity: (mm*2)/s 


$D1 = 4.1E-3 
/ Solid Diffusivity: (mm*2)/s 
$Ds = 4.1E-7 


/***#**e AKA TNITIAL ACTIVITY CONDITIONS ********* 
Peinitial Concentration: wt/ 


$C_tO = 10 
feiditial Activity: — 
$A_tO = 10 


/****K*KKKKK DERIVED QUANTITIES ****+*****+ 

/ Liquidus Temperature: C 

$Tl = ($m1*$C_t0)+$Tf 

/ Solidus Temperature: C 

$Ts = ($ms*$C_tO)+$T£ 

/***x***xK* INITIAL TEMPERATURE CONDITIONS ****«*4 «4 
/Initial Temperature of the Left Wall: C 
$T_x0_tO = $T1l 

/Initial Temperature of the Right Wall: C 
$T_xL_tO = $T1+15000 

/Temperature Gradient: C/mm 

$dTdx = ($T_xL_t0-$T_x0_t0)/$L 

/******k** PROBLEM DEFINED CONSTANTS ****#* 4 
/ Interface Velocity: mm/s 


$Vp = 0.01 

/ 

/*******K*K** PROBLEM DEFINITION ********** 
‘3 


PROBLEM (2-D , TRANSIENT , NONLINEAR , NOMOMENTUM , ENERGY , SPECIES , SPECIES=1) 
EXECUTION (NEWJOB) 
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PRINTOUT (NONE) 

DATAPRINT (NONE) 

/ 

/*** KKK CONTINUEM ENTITIES *&* kkk 4x 
/ 

ENTITY (FLUID, NAME="cusn" , PROPERTY="cusn" ) 
/ 

/***** KAKA BOUNDARY ENTITIES ****** kk KX 
/ 

ENTITY (PLOT , NAME="wall") 

ENTITY (PLOT , NAME="left") 

ENTITY (PLOT , NAME="right") 


/ 

/*** KKK SOLUTION STRATEGY ********** 

k 

SOLUTION (S.S.=20, VELCONV=0.001 , RESCONV=0.001) 
RELAXATION () 


mmo.) 0, 0.5; 0, 0, 0, 0.95 
TIME (BACK , DT=1E-6 , NSTEPS=1E5 , VARIABLE=1E-3 , INCMAX=1.5,DTMAX=0.1,TEND=1000) 
POSTPROC (NBLO=1) 
eet, £00 
/ 
/**x***K*K*KKK MATERIAL PROPERTIES *#* 4x 
/ 
DENSITY (SET="cusn" , CONSTANT=$rho) 
CONDUCTIVITY (SET="cusn" , CONSTANT=$k) 
SPECIFICHEAT (SET="cusn" , SUBROUTINE=7 ) 
mrex02t0, $1-xL_to, J eee 6 MO ie $cp, SEES ool 
DIFFUSIVITY (SET="cusn" , SUBROUTINE=6) 
$ml, $Tf, $Ds, $D1, $cap_s, $cap_l 
CAPACITY (SET="cusn" , SUBROUTINE=4) 
$ml, $Tf£, $cap_s, $cap_l 
/ 
/*****KKAKKX SOURCE CONDITIONS ********** 
/ 
SOURCE (SPECIES=1 , ENTITY="cusn" , SUBROUTINE=6) 
$ml, $Tf, $cap_s, $cap_l, $rho, $A_tO 
/ 
/**k**KKKKAK FLUID INITIAL CONDITIONS *****e%*xx 
/ 
ICNODE (TEMPERATURE, POLYNOMIAL=1 , ENTITY="cusn") 
Biexo- tO. $dldx,1,0,0 
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ICNODE(SPECIES=1 , CONSTANT=$A_t0O, ENTITY="cusn") 

/ 

/ KK TOP/BOTTOM WALL BOUNDARY CONDITIONS #4 KX 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="wall") 

/ Activity Flux = 0 

/ 

/x#x«eKKKA RIGHT WALL BOUNDARY CONDITIONS **##4KHHK# 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="right") 

/ Activity Flux = 0 

/ 

/**x*x*x*e KA LEFT WALL BOUNDARY CONDITIONS ****k* x 
A 

BCNODE (TEMPERATURE , FSUBROUTINE, ENTITY="left") 

/ Activity Flux = 0 


/ 

RENUMBER (PROFILE) 

END() 

y 

/$$$$$$$S$$$ FISOLV H$$$$HSH$$SH$ 
/ 

CREATE(FISOLV) 

/ 
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B.4 Change in Interface Velocity Restart File 


/$$$$$$$$$$ DESCRIPTION $$$$$$H$$$ 
/ Intial-Steady_State-Final Problem 
/ Time-Varying Temperature BC 

/ No Flux Activity BC 

/ Interface Velocity = 0.01 mm/s 

/ 

PeILEC ) 

1-D DIRECTIONAL SOLIDIFICATION OF COPPER-TIN SOLUTION 
/ 

/$$$$$$$$$$ FI-GEN $$$$$S$$S$ 
FE-GEN(. ) 

/**e RRR KAA POINTS RRR KKK 
POINT (ADD , COORDINATE, X=0 , Y=-2) 
POINT (ADD , COORDINATE, X=30 , Y=-2) 
POINT (ADD , COORDINATE, X=0 , Y=2) 
POINT (ADD, COORDINATE, x=30 , Y=2) 
[RK KRKKKA LINES #2 RRR K 

POINT (SELECT , ID=1) 

POINT (SELECT , ID=2) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=3) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=1) 

POINT (SELECT, ID=3) 
CURVE (ADD , SEGMENT ) 

POINT (SELECT , ID=2) 

POINT (SELECT , ID=4) 
CURVE (ADD , SEGMENT ) 

[**KKKKKKKKK EDGING *##**KKKKX 
CURVE (SELECT, ID=1) 

CURVE (SELECT , ID=2) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=150, RATIO=0 , 2RATIO=0 , PCENTER=0 ) 
CURVE (SELECT , ID=3) 
CURVE(SELECT , ID=4) 

MEDGE (ADD , SUCCESSIVE, INTERVAL=2 , RATIO=0 , 2RATIO=0 , PCENTER=0) 
[ **RRRKRKAK MESH #2 RRR 
CURVE(SELECT , ID=1) 
CURVE(SELECT , ID=4) 

CURVE (SELECT , ID=2) 
CURVE(SELECT , ID=3) 
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MFACE (WIRE, LABEL="cusn" , EDG1=1 , EDG2=1 , EDG3=1 , EDG4=1) 
MFACE(SELECT , ID=1) 

MFACE (MESH, MAP, ENTITY="cusn") 

/ 

END( ) 

y 

/$$$$$$$$$$ FI-BC $$H$$$$H$$H$$S 

/ 

BaI-BC( ») 

BGADD (SELECT , EDGE, ID=1) 

BGADD (SELECT , EDGE, ID=3) 

BGADD (ADD , EDGE, ENTITY="wall" , INCLUDE) 
BGADD (SELECT , EDGE, ID=2) 

BGADD (ADD , EDGE, ENTITY="right" , INCLUDE) 
BGADD (SELECT , EDGE, ID=4) 

BGADD (ADD, EDGE, ENTITY="left" , INCLDUE) 
END( ) 

/ 

/$$$$$$$$$$ FIPREP $$$$$SH$S$ 

/ 

FIPREP(  ) 

/ ****K*K*KKKKX MESH DIMENSIONS *&***& KKK 
/ Mesh Length: mm 


$L = 30 
/ Mesh Width: mm 
SH = 4 


/********** PHASE DIAGRAM CONSTANTS #44444 
/ Pure Copper Freezing Temperature: C 

$Tf = 1084.87 

/ Peritectic Temperature: C 

$Tp = 798 

/ Peritectic Alpha Concentration: C 

$Ca = 13.5 

/ Peritectic Beta Concentration: C 

$Cb = 25.5 

/ Liquidus Slope: C/wt% 

$ml = ($Tp-$Tf£)/$Cb 

/ Solidus Slope: C/wt% 

$ms = ($Tp-$Tf)/$Ca 

/ Segregation Coefficient: - 

$kO = $ml1/$ms 

[4 eKKKKK CAPACITY-ACTIVITY RELATIONSHIP ORR 
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/ Solid Capacity: - 


$cap_s = $k0O 
/ Liquid Capacity: - 
$cap_l = 1 


/********** MATERIAL CONSTANTS #44 4X 
/ Density: g/(mm73) 

$rho = 7.75E-3 

/ Thermal Conductivity: W/(mm*k) 

$k = 0.35 

/ Specific Heat : J/(g*K) 

$cp = 0.516 

/ Latent Heat of Fusion: J/g 

$Lf = 236.1290323 

/ Liquid Diffusivity: (mm*2)/s 


$D1 = 4.1E-3 
/ Solid Diffusivity: (mm*2)/s 
$Ds = 4.1E-7 


Jee INITIAL ACTIVITY CONDITIONS HOG 
/ Initial Concentration: wt/ 


$C_tO = 10 
Peinitial j\Activity: - 
$A_tO = 10 


/******KKAK DERIVED QUANTITIES *#** 4% 

/ Liquidus Temperature: C 

$T1l = ($m1*$C_t0O)+$Tf£ 

/ Solidus Temperature: C 

$Ts = ($ms*$C_t0O)+$T£ 

/********** TNITIAL TEMPERATURE CONDITIONS *#* ++ 
/Initial Temperature of the Left Wall: C 
$T_x0_tO = $Tl 

/Initial Temperature of the Right Wall: C 
$T_xL_tO = $T1+15000 

/Temperature Gradient: C/mm 

$dTdx = ($T_xL_t0-$T_x0_t0)/$L 

/********** PROBLEM DEFINED CONSTANTS ****#& «k++ 
/ Interface Velocity: mm/s 


$Vp = 0.02 

/ 

[| RRR KKK KKK PROBLEM DEFINITION 2K OK 2K OK 2K 2K OK 2K OK OK 
/ 


PROBLEM (2-D , TRANSIENT , NONLINEAR , NOMOMENTUM , ENERGY , SPECIES , SPECIES=1) 
EXECUTION (RESTART) 
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PRINTOUT (NONE) 
DATAPRINT (NONE) 
/ 
[RRR KK CONTINUEM ENTITIES **********x 
/ 
ENTITY (FLUID, NAME="cusn" , PROPERTY="cusn" ) 
/ 
/ ***** KKK BOUNDARY ENTITIES *********x 
/ 
ENTITY (PLOT , NAME="wall") 
ENTITY (PLOT , NAME="left") 
ENTITY (PLOT , NAME="right") 
‘ 
/***e KKK SOLUTION STRATEGY *********X 
/ 
SOLUTION (S.S.=20, VELCONV=0.001 ,RESCONV=0.001) 
RELAXATION () 
maemo GO. FO.5. 10, 0,0) -0.95 
TIME (BACK, DT=1E-6 , NSTEPS=1E5 , VARIABLE=1E-3 , INCMAX=1.5, DTMAX=0. 1 , TEND=2000) 
POSTPROC (NBLO=1) 
fetes. LOO 
i 
/***k*k*k*k*K*K*K* MATERIAL PROPERTIES ******k**x 
/ 
DENSITY (SET="cusn" , CONSTANT=$rho) 
CONDUCTIVITY (SET="cusn" , CONSTANT=$k) 
SPECIFICHEAT (SET="cusn" , SUBROUTINE=7) 
BemrOStOy olexl tO, oml, STf, Scp, SLi, $L 
DIFFUSIVITY (SET="cusn" , SUBROUTINE=6) 
$ml, $Tf£, $Ds, $D1, $cap_s, $cap_l 
CAPACITY (SET="cusn" , SUBROUTINE=4) 
$ml, $Tf, $cap_s, $cap_l 
/ 
/***k*KKKKAKA SOURCE CONDITIONS *********x 
/ 
SOURCE (SPECIES=1 , ENTITY="cusn" , SUBROUTINE=6 ) 
$ml, $Tf£, $cap_s, $cap_l, $rho, $A_tO 
/ 
/##e**k*** TOP/BOTTOM WALL BOUNDARY CONDITIONS *#*+4 4x 
/ 
BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="wall") 
/ Activity Flux = 0 
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/ 

/******KKKX RIGHT WALL BOUNDARY CONDITIONS *******k* 
/ 

BCNODE (TEMPERATURE , FSUBROUTINE, ENTITY="right") 

/ Activity Flux = 0 

/ 

/***** eK LEFT WALL BOUNDARY CONDITIONS ***xxxxxx* 
/ 

BCNODE (TEMPERATURE, FSUBROUTINE, ENTITY="left") 

/ Activity Flux = 0 


i 

RENUMBER (PROF ILE) 

END( ) 

/ 

/$$$$$$$$$$ FISOLV $$$$$$S$$$ 
/ 

CREATE (FISOLV) 

/ 
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Appendix C 


FIDAP User Subroutines 


C.1 USRBCN Subroutine 


SUBROUTINE USRBCN (VAL,NODE,IDF,TIME,SOL,NUMEQA,NDOF ,NUMNP,LDOFU, 


C 
C USER DEFINED BOUNDARY CONDITIONS 
C 
C VAL COMPUTED (SPECIFIED) BOUNDARY CONDITION 
C SOL GLOBAL SOLUTION VECTOR 
C NUMEQA = GLOBAL EQUATION NUMBER ARRAY 
C NDOF = ACTIVE NUMBER OF DEGREES OF FREEDOM 
C NUMNP = NUMBER OF NODAL POINTS 
C NODE NODE NUMBER OF B.C. 
C IDF DEGREE OF FREEDOM FOR NODE 
C LIME PisaliMe 
C LDOFU = ACTIVE DEGREE OF FREEDOM ARRAY 
C CONSTR = ARRAY OF SPECIFIED NONZERO BOUNDARY CONDITIONS 
C XYZ nodal coordinates 
C iflag flag for user to set (not equal to 0) if coordinates 
C are updated 
C NODEPR = reverse permutation array 
C node(external) = NODEPR(NODE) where NODE = internal node no. 
C 
#include "IMPLCT.COM" 
#include "PARUSR.COM" 
C 
PARAMETER(T£ = 1084.87D0) 
PARAMETER(Tp = 798.0D0) 
PARAMETER(Ca = 13.5D0) 
PARAMETER(Cb = 25.5D0) 


CONSTR, NODEPR, XYZ, iflag) 


PARAMETER(C_tO = 10.0D0) 
PARAMETER(dx = 5.0DO) 
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PARAMETER(Vp1 = 0.010D0) 
PARAMETER(Vp2 = 0.020D0) 


DIMENSION SOL(*) ,NUMEQA(NUMNP ,NDOF) , CONSTR(*) , LDOFU (*) 
DIMENSION NODEPR(*) ,XYZ(NUMNP, *) 
F_DOUBLEPRECISION GETSOL 
F_DOUBLEPRECISION GETSOLP 
ZRO = 0.DO 


slope_l = (Tp-Tf)/Cb 
slope_s = (Tp-Tf)/Ca 
seg = slope_l/slope_s 
Tl = (slope_1*C_t0)+Tf 
Ts = (slope_s*C_t0)+Tf 
Case 1: Initial Transient 
oex0_toO = [1 
tex tO -=.f1+125020D0 
Thermal gradient and freezing rate 
dhox 6 (Tex LO0-18x0_t0)/dx 
T_doti = Vp1*dTdx 
T_dot2 = Vp2*dTdx 


Freezing at Vpl 

VAGe=)(Cieximto-) x08G0) #XYZ (NODE, 1)/dx)+T.x0_t0-(T_dot1*TIME) 
Vertical Freezing at Vp1l 

VAL = ((T_xL_t0-T_x0_t0)+*XYZ(NODE, 2) /dx)+T_x0_t0-(T_dot1*TIME) 
Restart values, was at Vpi, now at Vp2 

VAL = ((T_xL_t0-T_x0_t0) *XYZ (NODE, 1) /dx)+T_x0_t0-(T_dot1*1000.0DO) 

& -(T_dot2*(TIME-1000.0D0) ) 


RETURN 
END 
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C.2 USRCAP Subroutine 


SUBROUTINE USRCAP (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP, 
1 DSDX,XYZL,PROP,TIME,NPTS,ndp,MNDP, IERR) 


USER DEFINED SPECIES CAPACITY 


NELT = GLOBAL ELEMENT NUMBER 

NE = LOCAL ELEMENT NUMBER 

NG = GROUP NUMBER 

COEF = SPECIES CAPACITY 

VARI = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS 

DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS 
NDFCD = 2(3) COORDINATES DIMENSION, ACCORDING TO THE DEFINITION IN 
LDOFU = pointer array for accessing vari and dvari information 
XYZL = X,Y,Z COORDINATES 

SHP = ELEMENT SHAPE FUNCTIONS 

DSDX = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION 

PROP = USER DEFINED PARAMETERS 

MNDP = FIRST DIMENSION OF SHAPE FUNCTION MATRICES 

TIME = TIME 

NPTS = NUMBER OF POINTS 


PROP(1) = Liquidus Slope (ml) 
PROP(2) = Pure Cu Freezing Temperature (Tf) 
PROP(3) = Solid Capacitance (cap_s) 


PROP (4) = Liquid Capacaitance (cap_1l) 
PROP(5) = Solidus Temperature (Ts) 


CoG) GGG) ei GG GG) eG) GGG) Ge Gm Gms Ge) a Ge) eG Ge eC Ce) eee 


#include "IMPLCT.COM" 
#include "PARUSR.COM" 
C 
C My Code Part 1 (Start) 
C Constants 
PARAMETER(maxdim = 100000) 
PARAMETER(NUMNODE = 4) 
PARAMETER(NUMPOINTS = 9) 
PARAMETER(NUMSIDES = 2) 
C My Code Part 1 (End) 
C FIDAP Defined Parameters 
DIMENSION COEF (NPTS) 
DIMENSION SHP(NPTS,MNDP) ,DSDX(NPTS,NDFCD,MNDP) , XYZL(NPTS , NDFCD) 
DIMENSION PROP(*) , VARI (NPTS,*) ,DVARI(NPTS ,NDFCD,*) , LDOFU(*) 
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C 


C My Code Part 2 (Start) 
C Matricies 


C 
C 
C 
C 
C 
C 
C 
C 10 
C 

C 

C 


C SHP 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


shp_1(NUMNODE , NUMNODE) 
temperature_node(NUMNODE) ,activity_node (NUMNODE) 
x_node (NUMNODE) , y_node (NUMNODE) 

T_CV (NUMNODE , NUMNODE) , A_CV (NUMNODE , NUMNODE) 
X_CV(NUMNODE , NUMNODE) , Y_CV (NUMNODE , NUMNODE) 
phi_CV(NUMNODE, NUMNODE) 
nodal_sum_of_on(NUMNODE) 

phi_O_CV (NUMSIDES , NUMNODE) 

volume_fraction (NUMNODE) 

initial_call (maxdim) 

previous (maxdim , NUMNODE) 


LOGICAL USRCAP_test / .FALSE. / 
SAVE initial_call, previous 


Initial Values for restarts 


if (initial_call(1).eq.0.and.NE.eq.1)then 
initial calhk@iyr = I 
open (67) 


do 10 


j = 1,750 


read(67,*) previous(j,1), previous(j,2), 
& previous(j,3), previous(j,4) 
continue 
close(67) 


endif 


Initial Calculations for pressure where NPTS = 1 


if (NPTS.eq.1)then 
if (initial_call (NE) .eq.0)then 
initial_call(NE) = 1 
COEF(1) = PROP (4) 


else 


COEF(1) = 0.25D0* (previous (NE, 1)+previous (NE, 2) 


endif 


RETURN 


endif 


& +previous(NE,3)+previous (NE, 4) ) 


Matrix Inverse 


A = SQRT(3.0D0)/2.0D0 
B = 1.0DO0+A 
C = 1.0D0-A 


shp_1(1,1) = B 
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Shing 1:(452) 
shp_1(1,3) 


shp_1(1,4) = 


shp_1(2,1) 
shp_1(2,2) 
shp_1(2,3) 
shp_1(2,4) 


shp_1(3,1) = 


shp_1(3, 2) 
shpsii3e3) 


shp_1(3,4) = 


shp_1(4,1) 
shp_1(4, 2) 


shp_1(4,4) 


= -0.5D0 
= -0.5D0 


= -0.5D0 


= -0.5D0 


= -0.5D0 
= -0.5D0 


=08 D0 


= -0.5D0 


C Find the solution at the nodes and their position 


30 
20 


do 20 n = 


1 , NUMNODE 


temperature_node(n) = 0.0DO 
activity_node(n) = 0.0DO 


x_node(n) = 0.0DO 
y_node(n) = 0.0DO 
ado 350-1 


= 1,NPTS 


temperature_node(n) = temperature_node(n)+(shp_1(n,i) 


& *VARI(i,LDOFU(KDT) ) ) 
activity_node(n) 
& *VARI(i,LDOFU(KDS+1) ) ) 


= activity_node(n)+(shp_1(n,i) 


x_node(n) = x_node(n)+(shp_1(n,i)*XYZL(i,1)) 


y_node(n) 


continue 


continue 


y_node(n)+(shp_1(n,i)*XYZL(i,2)) 


C Find solutions at the midpoints of the sides and at the center of the element 
CALL CORNERS (temperature_node,activity_node,x_node,y_node, 


eee CVA ASCV, X CV, Y_CV) 


C Determine which faces the interface passes through for each control volume 


do 40 i = 


1,NPTS 


GALES SIDES (TCV Gly 1) SAT CVG a) XeCV Clad) ayo Vly 
& PROP(1) ,PROP(2) ,NE,TIME, 
& phi_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i)) 
C Find liquid fraction for each control volume 
GALL AREA (X_CV(174) | YoCV.G 71); nodalesumsof2onG) ,phis0-CV(1,1)., 


& NE,TIME, 


& volume_fraction(i)) 
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C Find the capacity at each intergration point 
C Linear Interpolation 
COEF(i) = ((1.0D0-volume_fraction(i))*PROP(3) ) 
&+(volume_fraction(i)*PROP (4) ) 
C Save solution for pressure calculation 
previous(NE,i) = COEF(i) 
40 continue 
C File output for restart files 
C if (TIME.eq.44.2586D0) then 
c write(67,*) (previous (NE,i) ,i=1,4) 
C endif 
C My Code Part 2 (End) 
RETURN 
END 
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C.3. USRDIF Subroutine 


SUBROUTINE USRDIF (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP, 
1 DSDX,XYZL, PROP, TIME,NPTS ,ndp,MNDP, ICOEF , IERR) 


USER DEFINED SPECIES DIFFUSIVITY 


NELT = GLOBAL ELEMENT NUMBER 
NE = LOCAL ELEMENT NUMBER 
NG = GROUP NUMBER 
COEF = SPECIES DIFFUSIVITY followed by the 
TENSOR MULTIPLIER FOR NONISOTROPIC DIFFUSIVITY 
VARI = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS 
DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS 
NDFCD = 2(3) COORDINATES DIMENSION, ACCORDING TO THE DEFINITION IN 
LDOFU = pointer array for accessing vari and dvari information 
XYZL = X,Y,Z COORDINATES 
SHP = ELEMENT SHAPE FUNCTIONS 
DSDX = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION 
PROP = USER DEFINED PARAMETERS 
MNDP = FIRST DIMENSION OF SHAPE FUNCTION MATRICES 
TIME = TIME 
NPTS = NUMBER OF POINTS 


PROP(1) = Liquidus Slope (ml) 

PROP(2) = Pure Cu Freezing Temperature (Tf) 
PROP(3) = Solid Diffusivity (Ds) 

PROP(4) = Liquid Diffusivity (D1) 

PROP(5) = Solid Capacitance (cap_s) 

PROP(6) = Liquid Capacitance (cap_l) 
PROP(7) = Solidus Temperature (Ts) 


AZQsQeOtoeOsoOmo aan Ge OaG eG eG). GrGls GlGln Gla Gls Gls Gla Gla Glm GlmGo 


#include "IMPLCT.COM" 

#include "PARUSR.COM" 

C 

C My Code Part 1 (Start) 

C Constants 
PARAMETER(maxdim = 100000) 
PARAMETER(NUMNODE = 4) 
PARAMETER(NUMPOINTS = 9) 
PARAMETER(NUMSIDES = 2) 

C My Code Part 1 (End) 
DIMENSION COEF (ICOEF, 10) 
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DIMENSION 
DIMENSION 


SHP(NPTS ,MNDP) ,DSDX(NPTS ,NDFCD,MNDP) , XYZL(NPTS , NDFCD) 
PROP (*) , VARI (NPTS, *) ,DVARI(NPTS,NDFCD, *) , LDOFU(*) 


C My Code Part 2 (Start) 
C Matricies 


AIOVQt tara Qe Oro Qiame) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


shp_1(NUMNODE , NUMNODE) 
temperature_node(NUMNODE) , activity_node (NUMNODE) 
x_node (NUMNODE) , y_node (NUMNODE) 

T_CV (NUMNODE , NUMNODE) , A_CV (NUMNODE, NUMNODE) 
X_CV (NUMNODE , NUMNODE) , Y_CV (NUMNODE, NUMNODE) 
phi_CV (NUMNODE , NUMNODE) 
nodal_sum_of_on(NUMNODE) 
phi_O_CV(NUMSIDES , NUMNODE) 
volume_fraction(NUMNODE) 
initial_call(maxdim) 

previous (maxdim, NUMNODE) 


ENGICAE USRDIF_test / .FALSE.: / 
SAVE initial_call,previous 


Initial Values for restarts 


if (initial_call(1).eq.0.and.NE.eq.1)then 
initial_call(1) = 1 
open (68) 


do 10 


fF=919750 


read(68,*) previous(j,1), previous(j,2), 
& previous(j,3), previous(j,4) 
continue 
close(68) 


endif 


if (NPTS.eq.1)then 
if (initial_call (NE) .eq.0)then 
initial_call(NE) = 1 
COEF(1,1) = PROP(4) 


else 


COEF(1,1) = 0.25D0* (previous (NE, 1)+previous (NE, 2) 
& +previous(NE,3)+previous (NE, 4) ) 


endif 


RETURN 


endif 


C SHP Matrix Inverse 
A = SQRT(3.0D0)/2.0D0 
B = 1.0D0+A 
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C = 1.0D0-A 
shp_1(1,1) = B 
shp_i(i:,2) ="-0-5D0 
shp_1(1,3) = -0.5D0 
shp_1(1,4) =C 
shp_1(2,1) = -0.5D0 
shp_1(2,2) =C 
shpeie273)ie=88 
shp_1(2,4) = -0.5D0 
shp_1(3,1) = C 
shp.t(3}2)=*-0-5D0 
shp_1(3,3) = -0.5D0 
shp_1(3,4) = B 
shp_1(4,1) = -0.5D0 
shp_1(4,2) = B 
shp_1(4,3) =C 
shp_1(4,4) = -0.5D0 
C Find the solution at the nodes and their position 
do 20 n = 1,NUMNODE 
temperature_node(n) = 0.0DO 
activity_node(n) = 0.0DO 
x_node(n) = 0.0DO 
y_node(n) = 0.0DO 
do 30 i = 1,NPTS 
temperature_node(n) = temperature_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDT) )) 
activity_node(n) = activity_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDS+1) ) ) 
x_node(n) = x_node(n)+(shp_1(n,i)*XYZL(i,1)) 
y_node(n) = y_node(n)+(shp_1(n,i)*XYZL(i,2)) 
30 continue 
20 continue 
C Find solutions at the midpoints of the sides and at the center of the element 
CALL CORNERS (temperature_node,activity_node,x_node,y_node, 
eaneGV ALCV,;X CV, Y_-CV) 
C Determine which faces the interface passes through for each control volume 
do 40 i = 1,NPTS 
CARL OIDES CleCV Glia), AeGV Git) CV ie Ye CV Cd: 
& PROP(1),PROP(2) ,NE, TIME, 
& phi_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i)) 
C endif 
C Find liquid fraction for each control volume 
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CALL AREA(X_CV(1,i) ,Y_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i), 
& NE,TIME, 
& volume_fraction(i)) 
C Find the diffusivity at each integration point 
COEF (i,1) = ((1.0D0-volume_fraction(i) ) *PROP (3) ) 
&+(volume_fraction(i)*PROP (4) ) 
C Save solution for pressure calculation 
previous(NE,i) = COEF(i,1) 
40 continue 
C File output for restart files 
C if (TIME.eq.44.2586D0) then 
C write(68,*) (previous (NE,i) ,i=1,4) 
C endif 
C My Code Part 2 (End) 
RETURN 
END 
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C.4 USRSRPH Subroutine 


SUBROUTINE USRSPH (NELT,NE,NG,COEF,VARI,DVARI,NDFCD,LDOFU,SHP, 


QaoaQa@ae OuGuGie Ge Qua) QeG 1 GsGrG) CG), GlaG)=GiaGieQ 1Ou OOO g@ (ONG) IGl so eGlsG) Gia Gle Cl aG 


DSDX ,XYZL,PROP , TIME, NPTS ,ndp, MNDP, IERR) 


USER DEFINED SPECIFIC HEAT 


NELT = 
NE = 
NG = 
COEF = 
VARI = 
DVARI = 
LDOFU = 
XYZL = 
SH ae = 
DSDX = 
PROP = 
MNDP = 
TIME = 
NPTS = 


GLOBAL ELEMENT NUMBER 

LOCAL ELEMENT NUMBER 

GROUP NUMBER 

SPECIFIC HEAT 

ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS 
GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS 
pointer array for accessing vari and dvari information 
X,Y,Z COORDINATES 

ELEMENT SHAPE FUNCTIONS 

SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION 
USER DEFINED PARAMETERS 

FIRST DIMENSION OF SHAPE FUNCTION MATRICES 

TIME 

NUMBER OF POINTS 


Rubinstein Problem 


PROP (1) 
PROP (2) 
PROP (3) 
PROP (4) 
PROP (5) 
PROP (6) 
PROP (7) 
PROP (8) 


= Initial Temperature (T_inf) 

= Place Holder 

= Liquidus Slope (ml) 

= Pure Cu Freezing Temperature (Tf) 
= Specific Heat (cp) 

= Latent Heat of Fusion (Lf) 

= Space Holder 

= Solidus Temperature (Ts) 


Smith Problem 


PROP(1) = Left/Bottom Wall Temperature (T_x0_t0) 
PROP(2) = Right/Top Wall Temperature (T_xL_t0) 
PROP(3) = Liquidus Slope (ml) 
PROP(4) = Pure Cu Freezing Temperature (Tf) 
PROP(5) = Specific Heat (cp) 
PROP(6) = Latent Heat of Fusion (Lf) 
PROP(7) = Mesh Length (L) 
PROP(8) = Solidus Temperature (Ts) 

#include "IMPLCT.COM" 
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#include "PARUSR.COM" 
C 
C My Code Part 1 (Start) 
C Constants 
PARAMETER(maxdim = 100000) 
PARAMETER(NUMNODE = 4) 
PARAMETER(NUMPOINTS = 9) 
PARAMETER(NUMSIDES = 2) 
PARAMETER(delta_time_min = 1.0D-7) 
PARAMETER(denominator_min = 1.0D-9) 
C My Code Part 1 (End) 
DIMENSION COEF (NPTS) 
DIMENSION SHP(NPTS,MNDP) ,DSDX(NPTS,NDFCD,MNDP) ,XYZL(NPTS, NDFCD) 
DIMENSION PROP(*) , VARI(NPTS,*) ,DVARI(NPTS,NDFCD,*) , LDOFU(*) 
C My Code Part 2 (Start) 
C Matricies 
DIMENSION initial_call(maxdim) 
DIMENSION time_old(maxdim) ,time_new(maxdim) 
DIMENSION shp_i(NUMNODE, NUMNODE) 
DIMENSION temperature_node(NUMNODE) , activity_node (NUMNODE) 
DIMENSION x_node(NUMNODE) , y_node (NUMNODE) 
DIMENSION T_CV(NUMNODE, NUMNODE) , A_CV (NUMNODE , NUMNODE) 
DIMENSION X_CV(NUMNODE , NUMNODE) , Y_CV (NUMNODE , NUMNODE) 
DIMENSION phi_CV(NUMNODE, NUMNODE) 
DIMENSION nodal_sum_of_on(NUMNODE) 
DIMENSION phi_O_CV(NUMSIDES , NUMNODE) 
DIMENSION temperature_old_int (maxdim, NUMNODE) 
DIMENSION temperature_new_int (maxdim, NUMNODE) 
DIMENSION volume_fraction_old(maxdim, NUMNODE) 
DIMENSION volume_fraction_new(maxdim, NUMNODE) 
DIMENSION top(NUMNODE) , bottom (NUMNODE) 
DIMENSION previous (maxdim, NUMNODE) 
LOGICAL USRSPH_test / .FALSE. / 
SAVE initial_call,time_old,time_new 
SAVE temperature_old_int ,temperature_new_int 
SAVE volume_fraction_old,volume_fraction_new 
SAVE iteration_number,previous 


Initial Values for restarts 
if (initial_call(1).eq.0.and.NE.eq.1)then 
initrabecallid)s=e1 
open (63) 


PELL SULD IMS 
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open (21) 
open (66) 
do 10 j = 1,750 

time_new(j) = 44.2586D0 

read(63,*) temperature_new_int(j,1), 
& temperature_new_int(j,2), temperature_new_int(j,3), 
& temperature_new_int(j,4) 
read(21,*) volume_fraction_new(j,1), 
& volume_fraction_new(j,2), volume_fraction_new(j,3), 
volume_fraction_new(j,4) 

read(66,*) previous(j,1), previous(j,2), 

& previous(j,3), previous(j,4) 
10 continue 

close(63) 

close(21) 

close(66) 

iteration_number = 1 
endif 


Q_Q2OsO Oe Oar TOTO Orta Cer toe 
&> 


if (NPTS.eq.1)then 
if (initial_call (NE) .eq.0)then 
COEF (1) = PROP(5) 
else 
COEF(1) = 0.25D0* (previous (NE, 1)+previous (NE, 2) 
& +previous(NE,3)+previous (NE, 4) ) 
endif 
RETURN 
endif 
C 
C Check to see if it is the first time USRSPH is called for this element 
if (initial_call(NE) .eq.0)then 
initial_call(NE) = 1 
time_old(NE) = 0.0DO 
time_new(NE) = TIME 
dow20ti stL NPIS 
C Standard Prolems 
volume_fraction_old(NE,i) = 1.0DO 
volume_fraction_new(NE, i) 1.0DO 
C Rubinstein Problem Temperature Initial Conditions 
C temperature_old_int(NE,i) = PROP(1) 
C temperature_new_int(NE,i) = PROP(1) 
C Smith Problem Temperature Initial Conditions 
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C Horizontal Solidification 

C temperature_old_int (NE, i) 
C & /PROP(7))+PROP (1) 
C 
C 
C 


( (PROP (2)-PROP (1) ) *XYZL(i, 1) 


temperature_new_int (NE,i) ( (PROP (2)-PROP(1)) *XYZL(i, 1) 
& /PROP(7))+PROP (1) 
Vertical Solidification 
temperature_old_int(NE,i) = ((PROP(2)-PROP(1))*XYZL(i ,2) 
& /PROP(7))+PROP(1) 
temperature_new_int (NE, i) 
& /PROP(7))+PROP (1) 
20 continue 
iteration_number = 1 
endif 
C Check to see if this is a repeated time step 
if (TIME.1t.time_new(NE) ) then 
time_new(NE) = TIME 
endif 
C If at a new time step (non-iteration) update "old" values with "new" values 
if (TIME-time_new(NE) .gt.delta_time_min) then 
time_old(NE) = time_new(NE) 
time_new(NE) = TIME 
do 30 i = 1,NPTS 
temperature_old_int(NE,i) = temperature_new_int (NE, i) 
volume_fraction_old(NE, i) volume_fraction_new(NE,i) 
30 continue 
if (NE.eq.1)then 
iteration_number 
endif 
else 
if (NE.eq.1)then 
iteration_number = iteration_number+1l 
endif 
endif 
C SHP Matrix Inverse 


( (PROP (2)-PROP (1) )*XYZL(i , 2) 


I 
= 


A = SQRT(3.0D0)/2.0D0 
B = 1.0D0+A 

C = 1.0D0-A 
Sliped.cimet) = B 
shp_1(1,2) = -0.5D0 
shp_1(1,3) = -0.5D0 
shp_1(1,4) =C 
shpeic2) ty) t="-075D0 
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shp_1(2,2) =C 
shp_1(2,3) =B 
shp_1(2,4) = -0.5D0 
shp_1(3,1) =C 
shp_1(3,2) = -0.5D0 
shp_1(3,3) = -0.5D0 
shp_1(3,4) =B 
shp_1(4,1) = -0.5D0 
shp_1(4,2) =B 
shp_1(4,3) =C 
shp_1(4,4) = -0.5D0 
C Find the solution at the nodes and their position 
do 40 n = 1,NUMNODE 
temperature_node(n) = 0.0DO 
activity_node(n) = 0.0DO 
x_node(n) = 0.0DO 
y_node(n) = 0.0D0 
do 50 i = 1,NPTS 
temperature_node(n) = temperature_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDT) ) ) 
activity_node(n) = activity_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDS+1) )) 
x_node(n) = x_node(n)+(shp_1(n,i)*XYZL(i,1)) 
y_node(n) = y_node(n)+(shp_1(n,i)*XYZL(i,2)) 
50 continue 
40 continue 
C Find solutions at the midopoints of the sides and at the center of the element 
CALL CORNERS (temperature_node,activity_node,x_node,y_node, 
& T_CV,A_CV,X_CV,Y_CV) 
C Determine which faces the interface passes through for each control volume 
do 60 i = 1,NPTS 
SAELeoO LDR hGVCiet) pAGGVCin Xe CVl ea ny CV.Cl as 
& PROP(3) ,PROP(4) ,NE,TIME, 
feonie CV. 11), nodal sum-otson (1), phi OscV Cli) ) 
C endif 
C Find the liquid fraction for each control volume 
CAUDOAREA( X= CV (led) Y .CViClea se nodalasum=of_on(i).,phisQsCV (1.1), 
& NE,TIME, 
& volume_fraction_new(NE, i) ) 
C Update values of temperature at the integration points 
temperature_new_int(NE,i) = VARI(i,LDOFU(KDT) ) 
C Find the effective specific heat at each integration point 
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top(i) = (PROP(5)*(temperature_new_int (NE, i) 
&-temperature_old_int(NE,i)))+(PROP(6)*(volume_fraction_new(NE,i) 
&-volume_fraction_old(NE,i))) 

bottom(i) = temperature_new_int(NE,i)-temperature_old_int (NE, i) 

if (ABS (bottom(i)) .1t.denominator_min)then 

COEF (i) = PROP(5) 
else 
COEF (i) = ABS(top(i)/bottom(i)) 
endif 
C Save solution for the pressure solution 
previous(NE,i) = COEF(i) 
60 continue 
C File output for restart files 
C if (TIME.eq.44.2586D0) then 
C write(63,*) (temperature_new_int(NE,i) ,i=1,4) 
C write(21,*) (volume_fraction_new(NE,i) ,i=1,4) 
C write(66,*) (previous (NE,i) ,i=1,4) 
C endif 
C My Code Part 2 (End) 
RETURN 
END 
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C.5 USRSRC Subroutine 


SUBROUTINE USRSRC (NELT,NE,NG,SORCE,VARI,DVARI,NDFCD,LDOFU,SHP, 
1 DSDX ,XYZL,deter,PROP,TIME,NPTS,ndp,MNDP,IERR, 


IOPT) 
USER DEFINED SOURCE FOR ENERGY OR SPECIES EQUATIONS 


NELT = GLOBAL ELEMENT NUMBER 

NE = LOCAL ELEMENT NUMBER 

NG = GROUP NUMBER 

SORCE = HEAT OR SPECIES SOURCE (RETURNED VALUES) 


XYZL = X,Y,Z COORDINATES 


R*determinant 
SHP = ELEMENT SHAPE FUNCTIONS 


PROP = USER DEFINED PARAMETERS 

MNDP = FIRST DIMENSION OF SHAPE FUNCTION MATRICES 
TIME = TIME 

NPTS = NUMBER OF POINTS 

IOPT = 0 ENERGY EQUATION 

IOPT = N- TRANSPORT EQUATION FOR SPECIES N (0<N<16) 


PROP(1) = Liquidus Slope (ml) 

PROP(2) = Pure Cu Freezing Temperature (Tf) 
PROP(3) = Solid Capacity (cap_s) 
PROP(4) = Liquid Capacity (cap_l) 
PROP(5) = Reference Density (rho_0) 
PROP(6) = Initial Activitye (A_t0) 
PROP(7) = Velocity (Vp) 

PROP(8) = Diffusivity (D1) 

PROP(9) = Latent Heat of Fusion (Lf) 
PROP(10) = Initial Temperature (T_t0) 
PROP(11) = Solidus Temperature (Ts) 


2h (2) COP (OE LOY 2) Ley (Sy) eel Ce) Loe) tee) Leph (Pe er Tel (ey Teh (Pl (ei Kew eh (py (ey (ep) (OY Cor (eek (Py CP ley coh (Pye) 


#include "IMPLCT.COM" 
#include "PARUSR.COM" 
C 
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VARI = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS 
DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS 
LDOFU = pointer array for accessing vari and dvari information 


DETER = the value of the transformation determinant at 
each of the npts points. For axi-symmetric problem it is 


DSDX = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION 
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C My Code Part 1 (Start) 
C Constants 
PARAMETER (maxdim = 100000) 
PARAMETER (NUMNODE = 4) 

PARAMETER (NUMPOINTS = 9) 
PARAMETER(NUMSIDES = 2) 
PARAMETER(delta_time_min = 1.0D-7) 
C My Code Part 1 (End) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


SORCE(NPTS) 

SHP (NPTS ,MNDP) ,DSDX(NPTS ,NDFCD ,MNDP) , XYZL(NPTS , NDFCD) 
PROP (*) , VARI(NPTS,*) ,DVARI(NPTS,NDFCD, *) , LDOFU(*) 
DETER (NPTS) 


C My Code Part 2 (Start) 
C Matricies 


PLAT (PUMP IP) to) Corie 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


initial_call(maxdim) 

time_old(maxdim) ,time_new(maxdim) 
shp_1(NUMNODE, NUMNODE) 
temperature_node(NUMNODE) , activity_node (NUMNODE) 
x_node (NUMNODE) , y_node (NUMNODE) 
T_CV(NUMNODE, NUMNODE) , A_CV (NUMNODE , NUMNODE) 
X_CV(NUMNODE, NUMNODE) , Y_CV(NUMNODE, NUMNODE) 
phi_CV(NUMNODE, NUMNODE) 
nodal_sum_of_on(NUMNODE) 
phi_O_CV(NUMSIDES , NUMNODE) 
activity_old_int (maxdim, NUMNODE) 
activity_new_int (maxdim , NUMNODE) 
volume_fraction_old(maxdim , NUMNODE) 
volume_fraction_new(maxdim, NUMNODE) 
previous (maxdim , NUMNODE) 


EOGIGADSUSRSRG_test (/ /RALSE. %/ 

SAVE initial_call,time_old,time_new 

SAVE activity_old_int, activity_new_int 

SAVE volume_fraction_old,volume_fraction_new 
SAVE iteration_number,previous 


Initial Values for restarts 


if (initial_call(1).eq.0.and.NE.eq.1)then 
inibialscall ¢i)isml 
open (18) 
open (62) 


do 10 


j = 1,600 


time_new(j) = 1000.0D0 
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read(18,*) activity_new_int(j,1), activity_new_int(j,2), 
& activity_new_int(j,3), activity_new_int(j,4) 
read(62,*) volume_fraction_new(j,1), 
& volume_fraction_new(j,2), volume_fraction_new(j,3), 
& volume_fraction_new(j,4) 
10 continue 
close(18) 
close (62) 
iteration_number = 1 
endif 


if (NPTS.eq.1)then 
if (initial_call (NE) .eq.0)then 
SORCE(1) = 0.0DO 
else 
SORCE(1) = 0.25D0* (previous (NE, 1)+previous (NE, 2) 
& +previous(NE,3)+previous (NE, 4) ) 
endif 
RETURN 
endif 


Check to see if this is the first time USRSRC is called for this element 


if (initial_call(NE) .eq.0)then 
initial_call(NE) = 1 
time_old(NE) = 0.0DO 
time_new(NE) = TIME 
do 20 i = 1,NPTS 


Smith Initial Conditions 
activity_old_int(NE,i) = PROP(6) 
activity_new_int(NE,i) = PROP(6) 


volume_fraction_old(NE,i) = 1.0D0O 
volume_fraction_new(NE,i) = 1.0D0O 


Interface Initial Conditions 

activity_old_int(NE,i) = PROP(6)*(1.0D0+((1.0DO-PROP (3) ) 

& /PROP (3) *EXP(-PROP (7) *XYZL(i,1)/PROP(8) ))) 
activity_new_int(NE,i) = PROP(6)*(1.0D0+((1.0DO-PROP (3) ) 

& /PROP (3) *EXP (-PROP (7) *XYZL(i,1)/PROP(8) ) )) 
volume_fraction_old(NE,i) = 1.0DO 
volume_fraction_new(NE,i) = 1.0DO 

20 continue 
iteration_number = 1 
endif 
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C Check to see if this is a repeated time step 
if (TIME.1t.time_new(NE) ) then 
time_new(NE) = TIME 
endif 


C If at a new time step (non-iteration) update "old" values with "new" values 


if (TIME-time_new(NE) .gt.delta_time_min) then 
time_old(NE) = time_new(NE) 
time_new(NE) = TIME 
do 30 i = 1,NPTS 
activity_old_int(NE,i) = activity_new_int (NE, i) 
volume_fraction_old(NE,i) = volume_fraction_new(NE,i) 
30 continue 
if (NE.eq.1)then 
iteration_number = 1 
endif 
else 
if (NE.eq.1)then 
iteration_number = iteration_number+1 
endif 
endif 
C SHP Matrix Inverse 


A = SQRT(3.0D0)/2.0D0 
B = 1.0D0+A 
C = 1.0D0-A 


shp_1(1,1) = B 
shp_1(1,2) = -0.5D0 
shpei(193)e= =025D0 
shp_1(1,4) =C 
sbpel (21 )e=*-025D0 
shpst¢2 /2)i=06 
shppi(273)e= B 
shp_1(2,4) = -0.5D0 
shp_1(3,1) =C 
shp_1(3,2) = -0.5D0 
shp_1(3,3) = -0.5D0 
shp_1(3,4) = B 
shp_1(4,1) = -0.5D0 
shp_1(4,2) = B 
shp_1(4,3) =C 
shp_1(4,4) = -0.5D0 
C Find the solution at the nodes and their position 
do 40 n = 1,NUMNODE 
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temperature_node(n) = 0.0DO 
activity_node(n) = 0.0DO 
x_node(n) = 0.0DO 
y_node(n) = 0.0DO 
do 50 i = 1,NPTS 
temperature_node(n) = temperature_node(n)+(shp_i(n,i) 
& *VARI(i,LDOFU(KDT) ) ) 
activity_node(n) = activity_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDS+1))) 
x_node(n) = x_node(n)+(shp_1(n,i)*XYZL(i,1)) 
y_node(n) = y_node(n)+(shp_1(n,i)*XYZL(i,2)) 
50 continue 
40 continue 
C Find solutions at the midopoints of the sides and at the center of the element 
CALL CORNERS (temperature_node,activity_node,x_node,y_node, 
Perey, A» CV, X<CV., Y-CV) 
C Determine which faces the interface passes through for each control volume 
do 60 i = 1,NPTS 
CARUSOLDES CleCViC1 1), AP GV Cia enkecv Cle) ve cyl, 1), 
& PROP(1) ,PROP(2) ,NE, TIME, 
cepul cv Clad) nodal ssum_ofon(i) ,phivO.CV(14i)) 
C endif 
Find the liquid fraction for each control volume 
CALL AREA(X_CV(1,i),Y_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i), 
& NE, TIME, 
& volume_fraction_new(NE, i) ) 
C Update values of activity at the integration points 
activity_new_int(NE,i) = VARI(i,LDOFU(KDS+1) ) 
Find the source to add at each integration point 
NOTE: Fidap solves the equation rho*D(C)/Dt = rho*d(Diff*d(A) /dx)/dx 
Therefore, SORCE must be pre-multiplied by density 
SORCE(i) = -PROP(5)*(PROP(4)-PROP(3))*activity_old_int (NE, i) 
& *(volume_fraction_new(NE,i)-volume_fraction_old(NE,i)) 
& /(time_new(NE)-time_old(NE) ) 
C Save solution for pressure calculation 
previous(NE,i) = SORCE(i) 
60 continue 


Q 


Oar 


C File output for restart files 

C if (TIME.eq.1000.0D0) then 

S write(18,*) (activity_new_int(NE,i) ,i=1,4) 

C write(62,*) (volume_fraction_new(NE,i) ,i=1,4) 
C endif 
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C My Code Part 2 (End) 
RETURN 
END 
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C.6 USRVSC Subroutine 


SUBROUTINE USRVSC (NELT,NE,NG,VISC,VARI,DVARI ,NDFCD,LDOFU,SHP, 
i} DSDX,XYZL,PROP,TIME,NPTS ,ndp,MNDP, IERR) 


USER DEFINED VISCOSITY 


NELT = GLOBAL ELEMENT NUMBER 

NE = LOCAL ELEMENT NUMBER 

NG = GROUP NUMBER 

VISC = VISCOSITY 

VARI = ARRAY OF SOLUTION VARIBLES AT INTEGRATION POINTS 
DVARI = GRADIENTS OF SOLUTION VARIABLES AT INTEGRATION POINTS 
LDOFU = pointer array for accessing vari and dvari information 
XYZL = X,Y,Z COORDINATES 

SHP = ELEMENT SHAPE FUNCTIONS 

DSDX = SHAPE FUNCTION DERIVATIVES IN THE X,Y,Z DIRECTION 
PROP = USER DEFINED PARAMETERS 

MNDP = FIRST DIMENSION OF SHAPE FUNCTION MATRICES 

TIME .=,TIME 

NPTS = NUMBER OF POINTS 


PROP(1) = Liquidus Slope (ml) 

PROP(2) = Pure Cu Freezing Temperature (Tf) 
PROP(3) = Solid Viscosity (mu_s) 

PROP (4) = Liquid Viscosity (mu_1l) 

PROP(5) = Solidus Temperature (Ts) 


QiQ’@ sO GOtOs@ 1 #OTO tas t@ “Ar AsQr Oat OrO Oro rG 


#include "IMPLCT.COM" 
#include "PARUSR.COM" 
C 
C My Code Part 1 (Start) 
C Constants 
PARAMETER(maxdim = 100000) 
PARAMETER(NUMNODE = 4) 
PARAMETER(NUMPOINTS = 9) 
PARAMETER(NUMSIDES = 2) 
C My Code Part 1 (End) 
C FIDAP Defined Parameters 
DIMENSION VISC(NPTS) 
DIMENSION SHP(NPTS,MNDP) ,DSDX(NPTS,NDFCD,MNDP) ,XYZL(NPTS , NDFCD) 
DIMENSION PROP(*) , VARI(NPTS,*) ,DVARI(NPTS ,NDFCD, *) , LDOFU(*) 
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C My Code Part 2 (Start) 
C Matricies 


(@) (phe) Le) (al (@) Oh 121 Tel iek ley (>) 


DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


shp_1(NUMNODE , NUMNODE) 
temperature_node(NUMNODE) , activity_node (NUMNODE) 
x_node (NUMNODE) , y_node (NUMNODE) 

T_CV (NUMNODE, NUMNODE) , A_CV (NUMNODE, NUMNODE) 
X_CV(NUMNODE, NUMNODE) , Y_CV (NUMNODE, NUMNODE) 
phi_CV (NUMNODE , NUMNODE) 
nodal_sum_of_on(NUMNODE) 
phi_O_CV(NUMSIDES , NUMNODE) 
volume_fraction(NUMNODE) 
initial_call(maxdim) 

previous (maxdim , NUMNODE) 


LOGICAL USRVSC_test / .FALSE. / 
SAVE initial_call,previous 


Initial Values for restarts 


if (initial_call(1).eq.0.and.NE.eq.1)then 
initial_call(1) = 1 
open (69) 


do 10 


j = 1,750 


read(69,*) previous(j,1), previous(j,2), 
& previous(j,3), previous(j,4) 
continue 
close(69) 


endif 


if (NPTS.eq.1)then 
if (initial_call (NE) .eq.0)then 
initial_call(NE) = 1 
VISC(1) = PROP(4) 


else 


VISC(1) = 0.25D0* (previous (NE, 1)+previous (NE, 2) 
& +previous(NE,3)+previous (NE, 4) ) 


endif 


RETURN 


endif 


C SHP Matrix Inverse 


A = SQRT(3.0D0)/2.0D0 
B = 1.0DO+A 
C = 1.0D0-A 


enti Cle) =2B 
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shp_1(1,2) = -0.5D0 
shp_1(1,3) = -0.5D0 
shp_1(1,4) =C 
shp_1(2,1) = -0.5D0 
shp_1(2,2) =C 
shp_1(2,3) = B 
shp_1(2,4) = -0.5D0 
shoes, tyr=2¢ 
shp_1(3,2) = -0.5D0 
shp_1(3,3) = -0.5D0 
shp_1(3,4) = B 
shp_1(4,1) = -0.5D0 
shp_1(4,2) = B 
shp_1(4,3) =C 
shp_1(4,4) -0.5D0 
C Find the solution at the nodes and their position 
do 20 n = 1,NUMNODE 
temperature_node(n) = 0.0DO 
activity_node(n) = 0.0DO 
x_node(n) = 0.0DO 
y_node(n) = 0.0DO 
do 30 i = 1,NPTS 
temperature_node(n) = temperature_node(n)+(shp_1i(n,i) 
& *VARI(i,LDOFU(KDT) )) 
activity_node(n) = activity_node(n)+(shp_1(n,i) 
& *VARI(i,LDOFU(KDS+1) ) ) 
x_node(n) = x_node(n)+(shp_1(n,i)*XYZL(i,1)) 
y_node(n) = y_node(n)+(shp_i(n,i)*XYZL(i,2)) 
30 continue 
20 continue 
C Find solutions at the midpoints of the sides and at the center of the element 
CALL CORNERS (temperature_node,activity_node,x_node,y_node, 
elev, A CV xeCVYECV) 
C Determine which faces the interface passes through for each control volume 
do 40 i = 1,NPTS 
CALL SIDES (CI2CV Ci) FASGVidad) x eCV Cie) y2CV (1,1), 
& PROP(1) ,PROP(2) ,NE,TIME, 
& phi_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i)) 
C endif 
C Find liquid fraction for each control volume 
CALL AREA(X_CV(1,i),Y_CV(1,i) ,nodal_sum_of_on(i) ,phi_O_CV(1,i), 
& NE,TIME, 


106 


oe : 


‘vee Gal Syectd? caper eutee:4d oes ayver tsb bate 


¥n.2_ig@, (RA AG ee 





afe o@ bow netie @6) Dt Seales Se el Si 


gov. 7.06, 2V2 3,4 ee SP Pe 





ie rkei: mbt ete)’ (Capea © 4e7eeenS 
~ CLEVE els, a2) Cee nana, ¥.= fnjehen ge 


402 _#, SVE 





oa te O& ¢ 
»y (ay 0 ; 
—— apa ow RiR 
Toa C6 @edy; Tr 
en | vi 
. Som, (Bi8 oe 
ays te, 4 - 
we .o: © (die (aa 
sage dew cohwe af? ee ev trpien. ocr 
ee) eee «a 
mM & ¢ (el ebm lategee 
02.050 (e)e@en, 4710hiaa 
ow. = (anaes | 
Op 6 o-(eyabanug 
qra.et & BO 
svresegeed « (oiler equtetiqgnget © 
e( Craven, 2) 1 a 
si rrvitxe © (alabed, Qeavesne 7 
. (Us AL aay a 











eva lane | 
p)¥etijw ,ehna, ererawe goes) CRIREED) : 
evu_7 9 5:8 VOL 
sreayt es ips Ah 
1 


ow ny L8T 


ne 






















& volume_fraction(i) ) 
C Find the capacity at each intergration point 
VISC(i) = ((1.0D0-volume_fraction(i) )*PROP (3) ) 
&+(volume_fraction(i)*PROP (4) ) 
C Save solution for pressure calculation 
previous(NE,i) = VISC(i) 
40 continue 
C File output for restart files 
C if (TIME.eq.44.2586D0) then 
C write(69,*) (previous(NE,i) ,i=1,4) 
C endif 
C My Code Part 2 (End) 
RETURN 
END 
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